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Abstract 

Grain-scale  sorption  mass  transfer  is  an  important  process  that  must  be  considered 
when  predicting  clean-up  time  and  choosing  remediation  techniques  for  subsurface 
hazardous  waste  contamination.  Rate-limited  sorption  is  responsible  for  the  rebound 
effect,  where  remediated  groundwater  is  recontaminated  by  desorption.  Sorbed 
contaminants  are  not  available  for  microbial  degradation,  and  the  desorption  rate  may 
govern  the  effectiveness  of  natural  attenuation  by  biodegradation.  Grain-scale  sorption 
nonequilibrium  is  generally  attributed  to  diffusive  transport,  either  in  soil  organic  matter  or 
in  mineral  micropores.  Typically  used  sorption  mass  transfer  models  either  fail  to 
reproduce  long-term  slow  desorption  (first-order  models),  or  are  based  on  diffusion  in 
assumed  (often  spherical)  grain  geometries.  New  multisite  models  have  been  proposed 
that  incorporate  more  realistic  grain  geometries.  To  validate  these  models,  we  have 
conducted  sorption  rate  experiments  with  paraffin,  nylon,  and  porous  ceramic  spheres. 
These  synthetic  surrogate  soils  were  chosen  for  their  differing,  but  known,  sorption 
coefficients,  diffusion  coefficients,  and  geometries.  Experiments  were  conducted  in  batch 
systems  containing  only  a  single  material  and  size,  as  well  as  distributions  of  two  or  more 
materials  and  sizes.  We  tested  the  ability  of  the  model  to  simulate  the  behavior  of  these 
systems  and  to  fit  system  parameters  from  rate  data. 
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INVESTIGATION  OF  SORPTION  MASS  TRANSFER  MODELS 
USING  SYNTHETIC  SOILS 

1.  Introduction 

Overview  of  Research  Interest: 

Chemicals  released  into  the  subsurface  move  through  the  soil  and  often 
contaminate  the  groundwater.  Hydrophobic  organic  chemicals  (HOCs)  such  as 
chlorinated  solvents  and  fuel  components  released  at  many  Air  Force  installations  are  of 
concern  because  these  chemicals  may  possess  mutagenic,  teratogenic,  and/or  carcinogenic 
properties.  The  processes  that  affect  the  transport  and  fate  of  these  contaminants  must  be 
understood  in  order  to  assess  exposure  risks  to  humans  and  the  environment  and  to 
develop  efficient  and  cost-effective  remediation  strategies. 

The  transport,  distribution,  fate  and  biological  availability  of  many  HOCs  in  the 
groundwater  are  highly  dependent  on  sorption  by  aquifer  solids;  therefore,  understanding 
sorption  phenomena  between  solute  and  sorbent  is  critical.  Although  sorption  is  often 
considered  to  be  an  instantaneous  process  for  modeling  purposes,  sorption  and  desorption 
may  not  always  be  fast  compared  to  other  fate  or  transport  processes  such  as  advection 
and  dispersion.  Studies  indicate  the  time  the  aquifer  solids  are  exposed  to  the  solute  even 
at  low  flow  rates  may  not  be  long  enough  to  permit  complete  sorption  equilibrium 
(Roberts  et  al.,  1986;  Ball  and  Roberts,  1991a,  Brusseau  and  Rao,  1991a).  Equilibrium 
sorption  expressions  are  not  always  able  to  predict  the  results  observed  in  experimental 
studies  (van  Genuchten  et  al.,  1974).  Sorption  nonequilibrium  exists  in  part  because  mass 
transfer  between  moving  groundwater  and  soil  organic  matter  (SOM)  is  governed  by 
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diffusion.  This  physical  nonequilibrium  can  have  marked  effects  on  chemical  transport 
resulting  in  breakthrough  curves  (concentration  versus  time  profiles  at  a  point  in  space) 
that  have  earlier  arrivals  and  longer  “tails”  than  would  be  predicted  assuming  equilibrium 
sorption  (Weber  et  al.,  1991).  Nonequilibrium  sorption  within  individual  soil  grains  has 
been  found  to  greatly  affect  the  rate  of  biotransformation  (Ramaswami  et  al.,  1994)  by 
limiting  the  degree  to  which  the  contaminant  is  bioavailable  (Scow  and  Alexander,  1992). 
Bioavailibility  is  important  because  it  often  accounts  for  the  persistence  of  compounds 
that  are  biodegradable  and  might  otherwise  readily  decompose  (Alexander,  1994). 
Knowledge  of  sorption  phenomena  will  improve  the  ability  to  predict  contaminant  mobility 
and  bioavailibility.  By  improving  sorption  mass  transfer  models,  more  accurate  exposiure 
and  risk  estimates  for  groundwater  pollutants  and  improved  clean-up  operations  are 
possible.  New  models  developed  for  this  purpose  require  rigorous  testing  using  media  of 
known  properties  that  affect  sorption  sites  (Heyse,  1994).  Experiments  using  synthetic 
media  with  known  geometries  and  sorption  characteristics  must  be  made  to  test  model 
validity. 

Research  Objectives; 

The  purpose  of  this  experimental  study  is  to  test  a  theory  that  relates  sorption  rate 
to  geometry  of  sorbent  particles.  Specifically,  this  research  seeks  to  determine  if  a 
heterogeneous  distribution  of  diffusion  domains  can  be  described  by  one  average 
geometry,  and  to  determine  if  that  average  geometry  can  be  predicted  from  experimental 
sorption  rate  data. 
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n.  Review  of  the  Literature 


Sorption  Fundamentals: 

Sorption  is  the  process  by  which  the  solutes  in  the  liquid  phase  associate  with  the 
sorbent  or  solid  phase.  The  chemical  structure  of  a  molecule  and  the  nature  of  the  solid 
phase  play  large  roles  in  determining  the  extent  of  sorption.  Sorption  is  used  to  indicate 
both  adsorption  onto  a  two-dimensional  surface,  and  absorption  into  a  three-dimensional 
matrix  (Schwarzenbach  et  ah,  1993).  Adsorption  is  the  accumulation  occurring  at  an 
interface  while  absorption  is  the  partitioning  between  two  phases,  and  sorption  includes 
both  adsorption  and  absorption.  With  absorption  the  equilibrium  partition  coefficient  (Kp) 
quantifies  the  solute  distribution  in  the  aqueous  and  solid  phases. 

For  nonpolar  organic  contaminants,  sorption  is  dependent  on  the  hydrophobicity  of 
the  compound.  The  HOC  in  this  study  is  a  neutral  nonpolar  organic  compound  where  the 
partition  coefficient  is  represented  by  Equation  (1)  (Karickhoff,  et  al.,  1979): 

,  (1) 

with 

(2) 

P. 

where: 

X  =  sorbed  phase  concentration  in  the  organic  matter  and  the  micropores, 

Com  =  concentration  ofsorbate  associated  with  the  SOM, 
fom  =  weight  fraction  of  solid  which  is  SOM, 
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Cim  =  average  aqueous  immobile  phase  concentration, 

Qim  =  immobile  phase  porosity,  and 
Cl  =  concentration  of  solute  in  solution. 

For  hydrophobic  compounds,  the  fraction  of  organic  matter  of  the  subsurface  solid 
phase  is  the  dominant  soil  characteristie  affecting  sorption.  Studies  indicate  that  many 
HOCs  are  sorbed  primarily  by  organic  matter  (Kariekhoff  et  ah,  1979;  Chiou  et  ah,  1983; 
Gschwend  and  Wu,  1985;  Rutherford  et  ah,  1992).  Sorption  can  be  viewed  as 
partitioning  into  the  organic  matter  rather  than  as  sorption  onto  the  sorbent  surface 
(Kariekhoff,  1981;  Chiou  et  ah,  1979).  Chiou  et  ah  (1979,  1983)  have  described  SOM 
sorption  as  a  partitioning  process  due  to  the  polymer-like  stmeture  of  soil  organic  matter. 
Normalization  of  the  Kp  values  hy  foe  results  in  a  parameter,  organic  carbon  distribution 
coefficient  (Koc)  (Kariekhoff,  1981).  The  Koc  is  independent  of  the  soil  and  a  function 
only  of  the  chemical.  Kp  is  proportional  to  the/<,c  with  Koc  being  the  proportionality 
constant  as  shown  in  Equation  (3): 

(3) 

Joe 

Studies  to  support  the  relationship  shown  include  Miller  and  Weber  (1986)  who  noted 
signifieant  decreases  in  sorption  when  the  organic  content  of  the  soil  was  removed. 

Typical  values  for  foe  for  surface  soils  range  from  0.01  to  0.08  (Hamaker  and  Thompson, 
1972)  and  0.0002  to  0.01  for  alluvial  sand  aquifer  materials  (Schwarzenbach  and  Westall, 
1981). 
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Equation  (3)  demonstrates  that  Kp  can  be  computed  by  measuring  the/oc  of  the 
subsurface  material  and  selecting  a  value  of  Koc.  The  value  of  Koc  can  be  estimated  from 
the  octanol-water  partition  coefficient  Ka»  describes  the  partitioning  of  an  organic 
chemical  between  a  polar  phase  (water)  and  a  relatively  nonpolar  phase  (1-octanol). 
Karickhoff  et  al.  (1979)  proposed  that  the  octanol-water  partitioning  closely  parallels  the 
sorption  in  the  soil  system.  Since,  several  researchers  have  proposed  empirical 
expressions  to  relate  Koc  to  either  the  water  solubility  or  K^,  of  a  chemical  (Karickhoff  et 
al.,  1979;  Chiou  et  al.,  1983;  Schwarzenbach  and  Westall,  1981).  These  empirical 
relationships  assumed  the  porous  medium  to  be  saturated  with  water.  Deviations  in 
sorption  have  been  observed  when  the  soil  is  dehydrated  and  when  an  organic  solvent  is 
present  in  place  of  the  water.  Chiou  et  al.  (1985)  demonstrated  in  aqueous  systems  that 
the  water  molecules  compete  with  the  solute  for  sorption  sites  on  mineral  surfaces.  They 
found  that  as  the  medium  is  dehydrated,  more  of  the  organics  are  able  to  sorb  onto  the 
mineral  surfaces.  Also,  in  the  presence  of  organic  solvents,  the  solute  is  highly  soluble; 
thus,  the  sorption  of  the  solute  significantly  decreases. 
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Sorbent  Distribution; 


Sorption  of  solute  into  the  sorbent  matrix  is  not  fully  understood.  In  order  to 
predict  sorption  behavior,  it  is  important  to  conceptualize  the  mechanism  by  which  the 
sorbent,  the  soil  organic  matter,  retains  the  HOC  molecules.  SOM  has  been  described  to 
exist  in  three  geometries;  as  discrete  particles,  as  surface  coatings,  or  within  micropores  of 
mineral  grains  (Augustijn,  1993).  Each  of  the  three  organic  matter  geometries  is  likely  to 
exist  in  soils  with  the  most  abundant  controlling  the  sorptive  process.  Hydrophobic 
molecule  diffusion  into  SOM  has  been  used  to  explain  rate  limited  sorption.  Diffusive 
transport  is  caused  by  random  molecular  movement  of  solute  from  areas  of  high 
concentration  to  areas  of  low  concentration  (Crank,  1975).  Pick’s  first  law  quantifies 
diffusion  by  the  mathematical  equation  (Equation  (4)): 

J^-D-  (4) 

where: 

J  =  mass  flux  of  individual  components, 

D  =  diffusion  coefficient, 

C  =  concentration  of  individual  components,  and 
X  =  distance. 

The  magnitude  of  the  molecular  motion  is  determined  by  the  internal  energy  of  the 
solute  molecule,  dependent  on  temperature,  and  represented  by  the  diffusion  coefficient 
(Szecsody,  1988).  The  diffusion  of  an  aqueous  solute  in  the  subsurface  is  a  very  slow 
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process  with  diffusion  coefficients  on  the  order  of  10'^  or  10'®  cm^s  for  most  fuel 
constituents  (Schwarzenbach  et  al.,  1993). 

Rate  limited  sorption  may  be  explained  using  either  the  intraorganic  matter 
diffusion,  lOMD,  theory  (Karickhoff  and  Morris,  1985)  or  by  the  retarded  intraparticle 
diffusion  (RIPD)  theory  (MiUer  and  Weber,  1984, 1986;  Ball  and  Roberts,  1991b).  Each 
theory  was  developed  through  experimental  observations  under  differing  conditions  that 
might  serve  to  explain  the  differing  theoretical  conclusions. 

The  lOMD  theory  of  sorption  mass  transfer  likens  SOM  to  a  polymer  matrix.  In 
lOMD  nonequilibrium  is  caused  by  diffusion  limitations  into  the  polymer-like  matrix  of 
SOM.  Diffusion  path  lengths  in  the  SOM  matrix  can  be  high  enough  to  cause  significant 
mass  transfer  limitations.  Lee  et  al.  (1988)  investigated  the  movement  of  TCE  and  p- 
xylene  in  two  sand  aquifer  materials  and  concluded  that  the  nonequilibrium  was  indeed 
caused  by  diffusion  imitations  into  the  organic  matter  matrix  of  the  soils.  The  diffusion 
path  is  not  a  simple  straight  line  but  a  tortuous  path  between  the  polymer-like  humic 
molecules.  The  longer,  tortuous  path  results  in  longer  times  for  a  molecule  to  diffuse  into 
or  out  of  SOM.  Longer  path  lengths  would  indicate  increased  diffusion  limitations  and 
decreased  sorption  rates.  Brusseau  et  al.  (1991b)  and  Woodbum  et  al.  (1989)  proposed 
that  organic  cosolvents  increase  the  rate  of  sorption  due  to  their  ability  to  swell  SOM,  thus 
increasing  path  length.  The  researchers  conducted  batch  and  transport  experiments  to 
measure  sorption  of  polyaromatic  hydrocarbons  (PAHs)  on  Webster  soil  from  30/70  (v/v) 
and  50/50  (v/v)  methanol/water  mixtures.  The  studies  indicate  that  the  dominant 
mechanism  for  sorption  was  similar  to  the  mechanism  for  retention  of  PAHs  by  reverse 
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phase  liquid  chromatography  (RPLC).  RPLC  packing  material  is  known  to  swell  when 
used  as  a  sorbent  material,  and  when  used  with  a  cosolvent  mixture  the  rate  of  sorption  is 
increased.  Carroll  et  al.,  (1994)  investigated  diffusion  of  PCBs  from  the  swollen  and 
condensed  phases  of  the  Hudson  River  sediment  organic  matter  using  a  permeant/polymer 
diffusion  model.  Data  from  long-term  PCB  desorption  were  fit  and  equivalent  diameters 
for  swollen  and  condensed  humic  polymer  were  estimated.  A  3-fold  difference  in  diameter 
is  obtained  for  the  swollen  humic  polymer,  while  values  for  the  condensed  polymer  are 
virtually  identical.  The  larger  diameter  of  the  swollen  humic  polymer  results  in  longer  path 
lengths  which  would  indicate  increased  diffusion  limitations  and  decreased  sorption  rates. 

Gauthier  et  al.  (1987)  purport  the  major  cause  of  variations  in  the  observed 
coefficients  to  be  related  to  the  structural  characteristics  of  the  SOM.  They  found  the 
magnitude  of  the  Koc  values  correlated  strongly  with  the  degree  of  aromaticity  in  the 
humic  material.  Their  analysis  indicates  a  decrease  in  rate  for  polar  molecules  capable  of 
hydrogen  bonding  within  SOM  supports  lOMD. 

Retarded  intraparticle  diffusion  is  conceptually  described  as  sorption  in 
micropores.  Like  lOMD,  RIPD  has  also  been  postulated  as  a  possible  mechanism  for 
rate-limited  sorption  (Miller  and  Weber,  1984, 1986;  Wu  and  Gschwend,  1986;  Ball  and 
Roberts,  1991b).  In  contrast  to  the  flexible  pore  network  in  SOM  visualized  by  lOMD,  the 
micropores  of  RIPD  are  thought  to  be  rigid  and  permanent.  The  RIPD  geometry  is 
characterized  by  long,  tortuous  diffusion  path  lengths.  Slow  mass  transfer  is  the  result  of 
retardation  by  sorption  to  SOM  in  the  micropores  and  steric  hindrance  and  constrictivity 
between  the  advective  and  diffusive  domains  of  the  micropores.  Wu  and  Gschwend 
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(1986)  found  that  large  particles  showed  a  slower  approach  to  equilibrium  than  smaller 
but  otherwise  similar  particles  using  the  same  sorbate.  By  reducing  the  diffusive  path 
length  into  the  interior  of  the  particles  and  by  increasing  exposed  sorbent  surface  area, 
they  were  able  to  increase  sorption  rate.  Steinberg  et  al.  (1987)  mechanically  broke  soil 
particles  in  a  ball  mill  and  observed  an  increase  in  the  rate  of  desorption.  Apparently,  the 
reduction  in  pore  path  length  by  pulverization  induced  faster  rates.  Similarly,  Ball  and 
Roberts  (1991a)  analyzed  the  long-term  sorption  of  tetrachloroethene  (PCE)  and  1,2,4,5- 
tetrachlorobenzene  (TeCB)  by  aquifer  solids  and  found  increases  in  rate  with  particle 
pulverization  giving  evidence  of  RIPD. 

Both  the  RIPD  and  lOMD  theories  maintain  that  diffusion  can  be  used  to  explain 
rate  limited  sorption,  and  it  may  be  that  both  lOMD  and  RIPD  mechanisms  exist  in  the 
environment,  possibly  in  the  same  particle.  The  predominate  mechanism  in  soils  that  are 
high  in  SOM  appears  to  be  lOMD,  while  RIPD  may  predominate  in  soils  of  low  SOM 
content  or  soils  which  have  grain-scale  mineral  microporosity. 
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Sorption  Models; 


Mass  transfer  rate  limitations  are  often  modeled  using  either  diffusion  (second- 
order)  models  or  first-order  approximations  of  diffusion. 

First-Order  Models 

Early  physical  nonequilibrium  modeling  used  a  comparatively  simple  first-order 
rate  expression  to  describe  transfer  of  contaminant  between  regions  of  mobile  and 
immobile  water,  van  Genuchten  and  Wierenga  (1976)  described  the  mass  transfer  of 
solute  in  sorbing  porous  media  by  dividing  the  porous  media  into  four  domains:  (1)  mobile 
(turbulent)  water  located  in  the  larger  (inter-aggregate)  pores,  (2)  immobile  (stagnant) 
water  located  inside  aggregates  and  at  the  contact  points  of  aggregates  and/or  particles, 

(3)  dynamic  soil,  located  sufficiently  close  to  the  mobile  water  phase  so  that  equilibrium 
could  be  assumed  to  exist  between  solute  in  the  mobile  water  phase  and  the  soil  phase, 
and  (4)  stagnant  soil  where  equilibrium  could  be  assumed  to  exist  between  solute  in  the 
immobile  water  phase  and  the  soil  phase.  The  diffusive  rate-limitations  in  this  model 
comes  from  the  first-order  rate  expression  that  describes  transfer  of  contaminant  between 
mobile  and  immobile  water  regions. 

Based  on  these  assumptions.  Equation  (5)  was  developed  to  represent  transport 
between  the  domains  (van  Genuchten  arid  Wierenga,  1976) 

where: 

Cm  =  concentration  in  the  mobile  liquid  phase, 
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Cim  =  average  concentration  in  the  immobile  liquid  phase, 

0m  =  mobile  phase  porosity, 

0,m  =  immobile  phase  porosity, 

f  =  fraction  of  sorption  sites  that  are  in  contact  with  the  mobile  liquid  phase, 
pfe  =  bulk  soil  density, 

Kd  =  linear  sorption  distribution  coefficient. 

Dm  =  mobile  region  dispersion  coefficient, 

V  =  average  pore  water  velocity,  and 
t  =  time. 

The  mass  transfer  between  the  mobile  and  immobile  phases  was  considered  to  be  first- 
order  modeled  as  Equation  (6): 

_ 

where: 

a'  =  first-order  mass  transfer  coefficient. 

A  modification  of  the  approach  described  above  is  the  two-site  sorption  kinetic 
model  that  is  characterized  by  two  sorptive  sites.  Si  and  S2  (Brusseau  and  Rao,  1989). 
The  Si  sites  are  equilibrium  sites  where  the  solid  phase  rapidly  equilibrates  with  the 
solution  phase.  The  S2  sites  are  rate-limited  sites  where  sorption  between  the  aqueous  and 
solid  phase  is  described  by  a  first-order  rate  equation  (van  Genuchten  and  Wagenet, 

1989).  This  model  can  be  visualized  as 
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where: 

C  =  mass  of  solute  dissolved  in  the  liquid  phase 
S  =  mass  of  solute  sorbed  per  mass  of  sorbent 
F  =  fraction  of  sorbent  for  which  sorption  is  instantaneous 
k2  =  first-order  mass  transfer  rate  coefficient 
This  model  is  described  mathematically  by  Equations  (7)  -  (9): 

(7) 

Si=FKpC  (8) 

^  =  k,(K^C(}-F)-S,).  (9) 

Nkedi-Kizza  et  al.  (1989)  show  the  two-site  model  is  mathematically  equivalent  to 
the  van  Genuchten  and  Wierenga  (1976)  model.  The  three  parameters  needed  to  model 
two-site  mass  transfer  are  the  equihbrium  partition  coefficient,  Kp,  the  fraction  of 
equilibrium  sites,  F,  and  the  desorption  rate  coefficient  for  the  slow  sites,  ^2-  The  partition 
coefficient  can  be  obtained  independently  through  batch  isotherm  determinations  in  the 
laboratory.  F  and  are  obtained  by  fitting  models  to  experimental  data.  Although 
relatively  simply  in  representation,  first-order  rate  models  are  constrained  in  that  the 
partition  coefficient  must  be  derived  experimentally  or  estimated  before  rate  parameters 
can  be  fit;  also,  model  parameters  have  been  shown  to  be  dependent  on  pore  water 
velocity  (Brusseau  and  Rao,  1989;  Kookana  et  al.,  1993). 
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Diffusion  Models 

More  complex  diffusion  models  were  also  being  developed  to  describe  the  transfer 
of  solute  within  the  immobile  phase.  These  models  differed  from  the  ones  previously 
discussed  in  that  diffusion  was  considered  to  be  Fickian  rather  than  first-order  and  the 
immobile  region  concentrations  are  no  longer  assumed  to  be  uniform. 

Pick’s  second  law  of  diffusion  describing  the  liquid  phase  concentration 
distribution  in  the  immobile  zone  as  Equation  (10): 


dt 


1  a 

8^  38 


(D,8^ 


3C., 


im 


38 


(10) 


where: 

Cim  =  concentration  (at  path  length  8  and  time  t)  in  the  immobile  liquid  phase 
8  =  diffusion  path  length 
X  =  shape  factor 

Da  =  apparent  diffusion  coefficient. 

The  shape  is  defined  by  a  distribution  of  capacity  along  the  path  length  as  defined  by 
Equation  (11): 

/(8)  =  ^(8„^-l)^  ;{0<6  <5,3x1  (11) 

Generally  Pick’s  second  law  of  diffusion  is  shown  for  a  sphere  where  X  =  2  and  8max  =  the 
radius  of  the  sphere. 

Building  on  the  concepts  developed  by  van  Genuchten  and  Wierenga  (1976),  Rao 
et  al.  (1980a)  investigated  the  use  of  two-site  (bicontinuum)  diffusion  model  to  describe 
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nonequilibrium  breakthrough  curves.  They  compared  the  van  Genuchten  and  Wierenga 
(1976)  model  and  one  incorporating  advection,  dispersion,  and  diffusion  into  spherical 
aggregates.  While  both  models  gave  relatively  good  predictions  of  experimental  results, 
better  predictions  were  obtained  when  the  rate  of  mass  transfer  was  described  by  diffusion 
into  aggregates.  Experimental  observations  (Rao  et  al.,  1980a;  1982)  indicated  that  most 
sorption  data  exhibits  a  two-stage  approach  to  equilibrium,  with  a  rapid  initial  rate,  which 
accounts  for  approximately  20-50%  of  total  sorbed,  followed  by  a  much  slower  rate. 

Miller  (1984)  proposed  a  physical  nonequilibrium  model  that  incorporated  film 
transport  and  intraparticle  diffusion  as  the  source  of  the  nonequilibrium  breakthrough.  The 
model  by  Miller  (1984)  and  used  by  Miller  and  Weber  (1984)  was  developed  using  mass 
transfer  and  mass  balance  concepts.  The  dual  resistance  model  was  developed  to  describe 
sorption  as  a  series  of  mass  transfer  steps  involving  molecular  diffusion  through  a  film 
layer  outside  the  particle  followed  by  diffusion  into  the  particle  itself.  The  soil  phase 
concentration  of  solute  is  described  to  vary  with  time  as  a  function  of  the  radial  dimension. 

Crittenden  et  al.  (1986)  developed  a  Fickian  physical  nonequilibrium  model  similar 
to  that  of  Miller  (1984).  Their  model  was  based  on  the  presence  of  aggregates  in  the  soil 
that  caused  the  nonequilibrium  breakthrough  curves.  This  model,  the  dispersed  flow, 
pore,  and  surface  diffusion  model  (DFPSDM),  included  intraparticle  diffusion  both  in  the 
pore  space  and  along  the  pore  surfaces.  The  DFPSDM  included  constant  advective  flow, 
axial  dispersion  and  diffusion,  film  mass  transfer  resistance  from  the  mobile  to  the 
immobile  phase,  local  equilibrium  between  solute  adsorbed  onto  the  soil  matrix  and  solute 
in  the  intra-aggregate  stagnant  fluid,  and  surface  and  pore  diffusion  as  intra-stationary 
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phase  mass  transport  mechanisms.  As  with  previous  models,  the  soil  column  was 
idealized  as  containing  mobile  water  and  stationary,  uniform,  spherical  aggregates.  In 
addition  to  the  DFPSDM  model,  Crittenden  et  al.  (1986)  present  and  solve  several 
simplified  versions. 

Researchers  (Goltz,  1986;  Wu  andGschwend,  1986;  Ball  and  Roberts,  1991b; 
Young  and  Ball,  1995)  have  applied  Pick’s  second  law  of  diffusion  (Equation  (9))  where 
the  mechanism  for  diffusion  was  in  the  sorbed  phase  as  well  as  diffusion  in  the  water  filled 
pores.  The  models  assume  linear,  reversible  equilibrium  sorption  within  the  spherical 
immobile  geometry.  Although  these  models  allowed  good  predictive  capabilities  for  the 
laboratory  column  studies  relative  to  the  various  assumptions,  substantive  knowledge  of 
the  particle  geometry  is  necessary  to  invoke  the  models. 

New  Multisite  Models 

A  review  of  various  diffusion  and  first-order  approximations  of  the  diffusion 
models  to  date  indicates  that  existing  two-site  sorption  models  are  deficient  in  that  the 
model  parameters  are  dependent  on  pore  water  velocity  (Brusseau  and  Rao,  1989; 
Kookana  et  al.,  1993),  the  models  fail  to  reproduce  long-term  slow  desorption  (Karickhoff 
and  Morris,  1985;  Connaughton  et  at.,  1993),  and  existing  models  based  on  diffusion 
assume  grain  geometries  in  the  absence  of  physical  characterizations.  Thus,  Heyse  (1994) 
proposed  two  new  multisite  models  that  incorporate  more  reahstic  grain  geometries.  The 
models  are  the  multiple  sites  in  series  (MSS)  and  multiple  sites  in  parallel  (MSP)  models. 
The  models  build  on  the  concept  introduced  by  Connaughton  et  al.  (1993)  that  sorption 
sites  are  a  continuum  of  compartments,  each  defined  by  a  first-order  mass  transfer  rate. 
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The  models  both  assume  that  sorption  can  be  described  as  a  partitioning  process  and  the 
rate  limiting  step  is  mass  transfer  by  diffusion  into  either  an  organic  phase  or  in 
intraparticle  micropores  (Chiou  et  al.,  1979, 1983;  Brusseau  and  Rao,  1989;  Ball  and 
Roberts,  1991b,  Young  and  Ball,  1995).  Where  the  two-site  model  uses  an  average  rate 
constant  to  define  the  sizes  of  sorption  site  compartments,  the  MSP  uses  a  frequency 
distribution  for  the  mass  transfer  rates  as  was  used  by  Connaughton  et  al.  (1993).  The 
MSP  model  is  based  on  the  assumption  that  the  diffusion  path  between  the  bulk  aqueous 
phase  and  each  sorption  site  is  independent  of  the  diffusion  path  to  all  other  sites.  The 
advantages  of  the  MSP  model  are  that  it  is  relatively  simple  computationally,  it  allows 
variations  in  SOM  geometry  by  invoking  a  stochastic  approach,  and  the  number  of  fitting 
parameters  is  limited.  The  disadvantage  in  the  MSP  model  lies  in  the  assumption  that  all 
concentration  gradients  are  linear. 

The  MSS  model  describes  diffusion  through  a  sequence  of  compartments  using  a 
frequency  of  sorption  capacity  along  a  diffusion  path  (5)  to  define  the  size  of  the  sorption 
compartments.  The  advantage  of  the  MSS  model  is  that  while  the  concentration  gradient 
between  each  adjacent  compartment  is  linear,  the  concentration  gradient  over  the  entire 
diffusion  path  length  can  be  nonlinear  (e.g.  finite  difference  approach  to  diffusion).  Like 
the  MSP,  the  MSS  model  requires  only  a  limited  number  of  fitting  parameters.  The 
disadvantage  of  the  MSS  model  is  it  is  computationally  complex  and  an  average  particle  is 
assumed  -  although  the  geometry  is  not  limited  to  a  sphere.  This  research  effort  seeks  to 
investigate  the  validity  of  the  MSS  model,  while  the  MSP  model  will  be  investigated  under 
future  efforts. 
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The  MSS  model  should  better  describe  the  diffusion  process  into  uniformly  sized 
films,  spheres,  or  micropores  by  improved  accounting  of  the  abundance  of  sites  at  the 
soil/water  interface  and  deep  sites  in  the  soil  particle.  Based  on  the  approximately  shape 
of  soil  grains,  many  models  have  assumed  a  spherical,  uniformly  sized  aggregate  structure 
(van  Genuchten  and  Wierenga,  1976;  Rao  et  ah,  1980a  and  b;  Nkedi-Kizza  et  al.,  1989; 
Crittenden  et  al,  1986;  Parker  and  Valocchi,  1986;  Ball  and  Roberts,  1991b;  Rounds  et 
al.,  1993;  Young  and  Ball,  1995).  The  spherical  distribution  may  not  adequately  account 
for  the  abundance  of  sorption  sites  at  the  water/sorbent  interface  nor  the  existence  of 
sorption  sites  at  long  path  lengths.  The  spherical  model  is  inaccurate  in  part  because  real 
soil  consists  of  a  variety  of  distributions  (e.g.  films  layers,  cylinders,  spheres  and 
micropores).  The  effectiveness  of  the  MSS  model  lies  in  its  potential  abihty  to  better 
depict  real  average  soil  distribution  of  sorption  sites.  In  real  soil  the  distribution  of 
sorption  sites  is  unknown  because  the  soil  characteristics  (shape  factor  and  diffusion  path 
length)  and  the  effective  diffusion  coefficients  are  unknown.  The  MSS  model  overcomes 
these  unknowns  by  using  a  soil  sorption  characteristic  that  is  readily  assumed,  and  fitting 
the  remaining  soil  parameters.  For  example,  the  free  liquid  diffusivity  of  the  solute  in  the 
solvent  of  interest  could  be  used  in  place  of  the  unknown,  variable  effective  diffusion 
coefficient.  The  model  would  use  the  free  liquid  diffusivity  and  fit  the  effective  shape 
factor  and  diffusion  path  length  for  an  average  “virtual”  soil  particle. 

A  conceptual  MSS  model  is  depicted  along  with  the  conventional  two- site  model 
in  Figure  (1). 
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The  two-site  model  was  described  by  Equations  (7)  -  (9).  The  governing  equation 


for  a  batch  system  is  given  by  Equation  (12): 


dC  dS  ^ 


(12) 


where: 

Vi  =  volume  of  the  liquid,  and 
Ms  =  mass  of  the  solid. 

The  concentration  in  the  sorbed  phase  is  normalized  to  the  total  mass  of  sorbent. 
The  rate  of  change  of  sorbed  concentration  in  each  compartment  of  the  MSS  model  is 
described  by  Equations  (13a)  -  Equations  (13c).  If  j  =  1: 
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If  1  <j<N: 
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If  J  =  N: 
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The  fraction  of  sorption  sites  (Fj)  is  defined  by  frequency  distribution /f5j  (Equation  (14)): 

/^,=p7(8)rf8  (14) 

H 

where  the  frequency  distribution  of  diffusion  path  lengths  is  defined  by  Equation  (1 1). 

The  size  of  all  compartments  sums  to  unity  (Equation  (15)): 

(15, 

y=i 

and  the  sorbed  concentration  (S)  is  the  sum  of  the  sorbed  concentration  in  each 
compartment  (Equation  (16)): 

(16) 

7=1 

An  advantage  of  the  MSS  model  is  the  ability  to  describe  diffusion  in  nonspherical 
shapes  by  fitting  the  X  shape  factor  parameter.  This  distribution  can  describe  uniform 
coatings  (X=0),  cylinders  (k=l),  spheres  (^=2),  higher  order  geometries  (^  >  2),  as  well 
as  any  combinations  of  geometries  (X=any  real  number  greater  than  zero).  Distributions 
are  defined  by  only  two  parameters  {hmax  and  A,). 

Batch  sorption  experiments  were  used  to  generate  data  to  determine  whether  the 
MSS  sorption  mass  transfer  model  can  correctly  predict  the  known  geometry  of  the  nylon 
balls,  paraffin  balls  and  porous  ceramic  spheres  sorbent  material.  An  analytical  solution  to 
the  diffusion  into  various  geometries  in  a  batch  system  as  a  function  of  time  (Crank,  1975) 
was  used  to  confirm  the  accuracy  of  the  numerical  solution  from  the  MSS  model. 
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Sorption  Studies  Involving  Synthetic  Media: 


Natural  soils  vary  in  type  and  composition  of  mineral  and  organic  fraction  making 
them  less  desirable  for  use  in  validating  mass  transfer  models.  Researchers  have  used 
synthetic  media  of  known  and  controlled  compositions  to  analyze  mass  transfer  under 
conditions  of  greater  simplicity  than  investigations  in  soil  of  unknown  characteristics. 
Selection  of  a  surrogate  for  SOM  is  challenging  because  of  the  heterogeneous  nature  of 
SOM  in  situ,  and  the  choice  of  media  depends  on  the  sorption  process  of  interest; 
investigations  in  pore  diffusion  or  polymer-like  diffusion. 

Although  the  diffusion  of  molecules  in  the  polymer-like  SOM  is  still  under 
investigation,  much  is  known  about  the  diffusion  in  polymers  (Crank  and  Park,  1968). 
Therefore,  researchers  have  used  synthetic  organic  polymer  surrogates  in  investigating 
rate  limiting  behavior  with  diffusion  models.  Woodburn  et  al.  (1989)  used  bonded  phase 
packing  materials  of  chain  length  C-2-C-8,  which  is  commonly  used  in  reverse-phase 
liquid  chromatography,  to  investigate  the  thermodynamics  and  mechanisms  of  HOC 
sorption  to  these  materials.  Other  researchers  have  used  polyacrylamide  gel-exclusion 
chromatography  beads  of  saturated  sizes  of  75  to  150  micrometers  in  diameter  to 
determine  the  effects  of  diffusion  of  ’‘^C-labeled  phenol,  /?-nitrophenol  and  glutamic  acid 
on  the  kinetics  of  biodegradation  (Scow  and  Alexander,  1992;  Scow  and  Hutson,  1992). 
Rebhun  et  al.  (1992)  investigated  the  role  of  mineral  surface  adsorption  through  use  of  a 
synthetic  aggregate  coated  with  humic  acid  at  varying,  but  known,  amounts.  Carroll  et  al., 
(1994)  investigated  the  desorption  of  polychlorinated  biphenyls  (PCB)  from  Hudson  River 
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sediments  using  XAD-4,  a  high  BET  surface  area  polystyrene  resin  bead,  as  a  PCB 
adsorbent.  The  precursor  to  this  work  examined  spheres  of  nylon  66  of  diameter  0.238, 
0.318, 0.476, 0.635, 0.794  cm  and  similarly  sized  paraffin  beads  to  measure  mass  transfer 
parameters  of  the  sorptive  diffusion  of  HOCs  into  these  polymeric  materials  and  to 
develop  simulation  models  that  describe  their  behavior  (Heyse,  1994).  Recently  in 
analyzing  sorption  as  a  partitioning  process  combined  with  the  occurrence  of  adsorption, 
Xing  et  al.  (1996)  used  mbbery  hydrophobic  and  hydrophilic  polymers  and  a  glassy 
polymer  as  model  sorbent  materials.  They  also  used  a  mesoporous  amorphous  silica  as  a 
model  for  hydroxylated  mineral  silica. 

Likewise  researchers  have  used  surrogate  media  to  investigate  pore  diffusion  as  a 
cause  of  rate  limiting  mass  transfer.  Rao  et  al.  (1980b)  used  porous  ceramic  spheres  made 
of  fired  kaolinite  clay,  an  aluminosilicate  mineral  of  formula  Al2(0H)4Si205,  of  size  0.55 
and  0.75  cm  in  radius.  They  measured  mass  transfer  parameters  for  the  nonadsorbed 
solutes  and  and  used  the  parameters  in  the  development  of  diffusion  simulation 
models.  Scow  and  Alexander  (1992)  and  Scow  and  Hutson  (1992)  used  the  porous 
spheres  courtesy  of  Rao  et  al.  (1980b)  in  addition  to  the  polyacrylamide  gel-exclusion 
chromatography  beads  discussed  above  to  determine  the  effects  of  pore  diffusion  of  ^'Re¬ 
labeled  phenol,  /7-nitrophenol  and  glutamic  acid  on  the  kinetics  of  biodegradation.  In  a 
study  to  investigate  the  effects  of  pore  size,  particle  size  and  lack  of  microporosity,  Farrell 
and  Reinhard  (1994)  used  silica  gels,  glass  beads  and  clay  montmorillonite  to  analyze  TCE 
desorption.  The  differing  media  allowed  targeted  investigations  into  desorption 
phenomena  of  meso- ,  micro-  and  nonporous  solids. 
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The  use  of  synthetic  media  gives  researchers  more  control  to  investigate  sorption 
properties  on  media  of  known  characteristics  including  diffusion  coefficients,  distribution 
characteristics  and  geometry.  Data  gained  from  studies  on  synthetic  media  improve  the 
researcher's  ability  to  interpret  investigations  with  real  soil. 
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ni.  Methodology 


Materials: 

Solid  Phases 

The  sorbent  materials  in  this  study  were  chosen  to  allow  investigation  of  diffusive 
transport.  The  synthetic  polymers  were  chosen  for  their  ability  to  substitute  for  modeling 
mass  transfer  rate  limitations  of  SOM  and  for  their  differing,  but  known,  sorption 
coefficients,  diffusion  coefficients,  and  geometries.  The  materials  include  paraffin  wax 
(Aldrich  Chemical  Co.,  CAS  registry  #  8002-74-2),  nylon  66  spheres  of  sizes  0.238, 
0.318,  0.476,  0.635,  0.794  cm  diameter  (Small  Parts  Inc.),  and  porous  ceramic  spheres 
made  of  fired  kaolinite  clay,  an  aluminosilicate  mineral  of  formula  Al2(0H)4Si205,  of  sizes 
0.55  and  0.75  cm  in  radius  (Rao  et  al.,  1980b). 

Paraffin  wax  spheres  were  prepared  from  paraffin  chunks  as  follows.  A  small 
quantity  of  paraffin  was  placed  in  a  beaker  and  heated  just  to  melting.  The  liquid  paraffin 
was  then  poured  into  a  small  plastic  mold  which  upon  solidification  and  removal  formed 
small  paraffin  plugs.  The  plugs  were  then  formed  into  spheres  by  pressing  and  turning  a 
paraffin  plug  through  a  0.483  cm  diameter  hole  in  a  stainless  steel  template.  This  was 
completed  multiple  times  until  a  uniform  paraffin  sphere  of  0.476  cm  diameter  was 
formed. 

Solute 

Anthracene  (99%,  Aldrich  Chemical  Co.,  CAS  #  120-12-7),  a  component  of 
petroleum  fuels  and  coal  tar,  is  an  polyaromatic  hydrocarbon  of  molecular  weight  178.2 
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g/mol.  Anthracene  has  good  adsorption  at  254  nm  and  is  easily  analyzed  with  HPLC 
ultraviolet  detection.  The  compound  is  stable,  nonpolar,  and  of  low  volatility  making  it 
ideal  for  laboratory  sorption  studies. 

Liquid  Phases 

The  liquid  phases  of  all  experiments  consisted  of  a  1:1  (v/v)  solution  of  methanol 
(99.9  +  %,  Aldrich  Chemical  Co.,  CAS  #  67-56-1)  and  distilled  water  (Aldrich  Chemical 
Co.,  CAS  #  7732-18-5)  with  differing  concentrations  of  anthracene  (<  0.9  mg/L).  The  1:1 
methanol  and  water  ratio  was  chosen  to  reduce  sorption  by  vial  glass  and  Teflon  covered 
septa  and  to  accelerate  the  sorptive  process  under  laboratory  conditions.  To  protect 
anthracene  from  photodegradation  aU  liquid  phases  were  stored  in  amber  glassware. 
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Experiments; 


Isotherm 

Nylon  was  prepared  by  independently  crushing  nylon  spheres  of  diameter  0.794 
and  0.318  cm  by  first  placing  them  in  a  flask  of  liquid  nitrogen.  Once  frozen,  the  spheres 
were  removed  and  placed  between  a  durable  cloth  and  crushed  with  a  hammer  on  a  hard 
surface.  Four  initial  concentrations  of  the  liquid  phase  were  used;  0.81, 0.61, 0.41,  and 
0.20  mg/L  of  anthracene.  For  each  size  of  nylon  sphere,  approximately  0.07  grams  of  its 
crushed  form  was  weighed  and  placed  into  4  mL  amber  vials.  Each  concentration  of  stock 
solution  was  added  to  each  type  of  crushed  nylon  sphere,  0.794  and  0.318  cm,  with  one 
duplicate  each.  In  addition  four  vials  were  filled  with  each  concentration  of  stock  solution 
for  the  quantitative  determination  of  the  initial  concentration  of  these  solutions.  Two  vials 
containing  0.07  grams  of  crushed  nylon  of  each  size  were  filled  with  analyte  free  1:1  liquid 
phase  to  act  as  blanks.  Isotherm  experiment  set  up  data  is  included  in  Appendix  B. 

The  vials  were  agitated  continuously  on  a  shaker  table  to  ensure  that  bulk  liquid  was 
completely  mixed.  The  solutions  have  not  reached  equilibrium  to  date. 

Rate  Studies 

The  rate  of  sorption  was  measured  using  different  types,  sizes,  and  combinations 
of  the  solid  phases.  All  experiments  were  performed  in  batch  by  placing  the  solid  phase 
and  anthracene  containing  liquid  phase  in  Teflon  capped  amber  vials  or  bottles.  The 
containers  were  agitated  and  sampled  periodically  to  monitor  the  change  in  liquid  phase 
anthracene  concentration  using  HPLC  as  the  batch  system  reached  equilibrium.  Blanks 
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consisting  of  anthracene  free  liquid  phase  and  solid  phase  were  used  to  ensure  that  there 
was  no  interference  associated  with  the  solid  phase. 

Nylon  66 

Batch  sorption  rates  were  measured  for  five  sizes  of  nylon  spheres;  0.238, 0.318, 
0.476, 0.635, 0.794  cm  diameter  spheres  as  well  as  for  a  mixture  of  all  five  sizes.  Nylon 
spheres  were  weighed  to  be  as  close  to  0.57  grams  as  possible  and  placed  in  40  mL  amber 
vials,  with  the  number  of  balls  weighed  noted.  Each  size  nylon  sphere  including  the  mixed 
sizes  was  run  in  triplicate.  To  each  vial,  approximately  40  mL  of  0.81  mg/L  anthracene 
liquid  phase  was  added  then  weighed  using  an  analytical  balance.  In  addition,  two  blanks 
were  made  consisting  of  mixed  sizes  of  nylon  spheres  with  analyte  free  1:1  liquid  phase. 
Three  4  mL  vials  were  filled  with  the  0.81  mg/L  anthracene  liquid  phase  to  determine 
quantitatively  the  concentration  of  the  initial  condition  and  to  act  as  a  laboratory  control. 

A  summary  of  nylon  66  experimental  conditions  is  included  in  Appendix  B. 

For  each  type  of  solid  phase,  two  of  the  three  vials  were  sampled  approximately 
every  ten  days  by  extracting  not  more  than  one  milliliter  of  the  liquid  phase  and  analyzing 
by  HPLC.  Vials  were  weighed  after  sampling  to  determine  mass  of  removed  liquid  phase. 
As  a  control,  the  third  vial  remained  unsampled  until  it  was  detemuned  equilibrium  was 
reached.  In  addition,  the  initial  condition  vials  were  sampled  repeatedly  over  the  course  of 
the  experiment  to  investigate  possible  loss  of  anthracene  mass  due  to  volatilization, 
degradation,  and/or  sorption  on  vial  glass  and  septa. 
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Porous  Ceramic  Spheres 

Batch  diffusion  rates  were  independently  measured  for  porous  ceramic  spheres  of 
0.55  and  0.75  cm  radii.  The  spheres  were  placed  in  a  vacuum  flask  with  anthracene  free 
liquid  phase  and  allowed  to  become  saturated  under  reduced  pressure  over  a  24  hour 
period.  Eight  0.55  cm  radii  sized  spheres  that  were  placed  in  a  40  mL  amber  vial  with  19 
mL  of  0.80  mg/L  anthracene  liquid  phase.  Six  0.75  cm  radii  sized  spheres  were  placed  in 
a  250  mL  amber  bottle  with  35  mL  of  0.80  mg/L  anthracene  liquid  phase.  Both  of  these 
systems  were  prepared  in  duplicate.  A  blank  consisting  of  the  same  conditions  as  the  0.55 
cm  sized  spheres  was  placed  in  another  40  mL  amber  vial  with  19  mL  of  anthracene  free 
hquid  phase.  Each  vial  and  bottle  were  capped  with  Teflon  lined  septa  and  gently  agitated 
by  hand  periodically.  Each  container  was  sampled  total  duration  of  about  48  hours  by 
extracting  about  1  mL  and  analyzing  by  HPLC.  Vials  were  weighed  after  sampling  to 
determine  mass  of  removed  hquid  phase.  A  summary  of  the  porous  ceramic  sphere  (SB 
and  LB,  small  and  large  beads  respectively)  experimental  conditions  is  included  in 
Appendix  B. 

Mixed  Porous  Ceramic  and  Paraffin  Spheres 

Batch  diffusion  rates  were  also  measured  for  a  mixture  of  media,  which  contained  both 
sizes  of  ceramic  spheres  and  paraffin  spheres  of  0.249  cm  radius.  A  duphcate  pair  of 
conditions  was  prepared  consisting  of  six  0.55  and  three  0.75  cm  radii  sized  saturated 
spheres  as  well  as  seven  paraffin  spheres  that  were  placed  in  40  mL  amber  vials  with  31 
mL  of  0.80  mg/L  anthracene  hquid  phase  each.  In  addition,  a  blank  consisting  of  the  same 
amount  and  type  of  solid  phase  was  placed  in  another  40  mL  amber  vial  but  anthracene 
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free  liquid  phase  was  added.  Each  vial  was  capped  with  Teflon  lined  septa  and  gendy 
agitated  by  hand  periodically.  Each  container  was  sampled  for  total  duration  of  four 
weeks  by  extracting  about  1  mL  and  analyzing  by  HPLC.  Vials  were  weighed  after 
sampUng  to  determine  mass  of  removed  liquid  phase.  A  summary  of  the  mixed  paraffin 
and  porous  ceramic  sphere  (PMB)  experimental  conditions  is  included  in  Appendix  B. 
Analytical  Methods 

Quantitative  determinations  were  accomplished  using  a  Hewlett  Packard  1090L, 
High  Performance  Liquid  Chromatograph  (HPLC).  The  operating  conditions  consisted  of 
a  100  microliter  sample  size,  100%  methanol  as  mobile  phase,  mobile  phase  velocity  of  0.4 
mL/min,  and  oven  temperature  of  32  degrees  Celsius.  Separations  were  accomplished 
using  a  150  X  4.6  mm  Alltech,  Adsorbosphere  analytical  column,  with  a  CIS  bonded 
phase,  and  5  micron  packing.  Detection  was  accomplished  using  a  filter  photometric 
detector  set  at  254  nm.  The  instrument  was  calibrated  using  a  series  of  standards  of 
different  concentration.  The  standard  stock  solution  was  made  by  weighing  a  known 
amount  of  anthracene  in  a  100  mL  class  A  volumetric  flask  and  diluting  to  the  mark  with 
methanol.  A  series  of  standards  of  varying  concentrations  was  made  by  diluting  the 
standard  stock  solution  with  methanol  using  class  A  pipets.  Instrumental  controls  were 
performed  prior  to  analyzing  experimental  samples  and  again  after  every  ten  experimental 
samples  to  verify  instmment  integrity. 
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IV.  Data  Analysis  and  Model  Simulation 


Experimental  Data  Analysis: 

Mass  balance  was  used  to  determine  sorbed  phase  concentrations  of  the  solid 
phase  for  each  experiment.  The  loss  of  anthracene  mass  due  to  volatilization,  degradation, 
and/or  sorption  on  vial  glass  and  septa  was  assumed  to  be  negligible  but  was  be  quantified 
using  data  obtained  from  the  blanks  and  controls. 

Nylon  66  Sphere  Rate  Study 

Because  each  vial  was  sampled  repeatedly  by  removing  a  small  aliquot  for  analysis, 
the  boundary  conditions  changed  following  each  sampling  event  To  account  for  the  loss 
of  mass  from  sampling,  the  following  relation  (Equation  (17))  was  used  in  determining  the 
solid  phase  concentration: 


)1  -Cl(n)Vl(n) 


(17) 


Experimental  data  was  plotted  as  Cj  /C/  versus  time  in  days  as  seen  in  Figure  (2).  When 
Cs  is  calculated  on  a  partitioning  basis  (Equation  (17)),  the  value  of  Cj  ICi  should 
eventually  approach  the  equilibrium  Kp  if  the  sorption  mechanism  is  partitioning.  For  a 
system  involving  the  same  solute  and  sorbent,  G  ICi  should  eventually  approach  the  same 
equilibrium  Kp.  At  120  days  into  the  experiment  it  is  obvious  that  either  the  media  is  not 
at  equilibrium,  or  the  sorption  mechanism  is  not  partitioning  if  it  is  at  or  near  equilibrium. 
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Figure  2.  Uptake  of  Anthracene  by  Nylon  66  (pardoning  basis). 
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Figure  3.  Uptake  of  Anthracene  by  Nylon  66  (adsorption  basis). 


Conversely,  when  Cs  is  calculated  per  surface  area  of  the  solid  phase,  plotting  the 
value  of  Cs  ICi  over  time  would  approach  an  adsorption  coefficient,  Kad,  if  the  sorption 
mechanism  is  adsorption.  Calculating  Cs  per  surface  area  solid  phase  (Equation  (18)): 

n 

= - - - - .  (18) 

The  uptake  data  is  presented  on  an  adsorption  basis  in  Figure  (3).  The  data  are  grouped 
closer  together  than  in  Figure  (2).  This  may  indicate  that  sorption  of  anthracene  by  nylon 
66  is  an  adsorption  process.  It  could  also  mean  that  sorption  of  anthracene  by  nylon  66  is 
a  very  slow  partitioning  process,  since  the  early  stages  of  partitioning  would  look  similar 
to  adsorption. 

Figure  (4)  represents  paraffin  rate  data  from  deVenoge  (1996)  that  is  plotted  as 
Cj/C/  versus  time  in  days  with  Cs  calculated  on  a  partitioning  basis  as  defined  in  Equation 
(17).  As  expected  for  partitioning  media,  the  value  of  C,  /Ci  for  the  paraffin  over  time 
approached  the  same  equilibrium  Kp  for  all  sphere  sizes.  When  the  paraffin  data  are 
plotted  on  an  adsorption  basis,  all  sphere  sizes  do  not  approach  the  same  equilibrium  Kad 
(Figure  (5)).  These  data  clearly  demonstrate  that  the  sorption  mechanism  for  anthracene 
by  paraffin  is  partitioning. 

An  argument  in  favor  of  sorption  of  anthracene  by  nylon  66  by  partitioning  can  be  made 
by  examining  the  apparent  rate  at  which  the  anthracene  is  sorbed  by  sphere  size.  If  the 
sorption  mechanism  was  on  an  adsorption  basis,  the  larger  spheres  would  have  the  fastest 
uptake.  Figure  (3)  clearly  shows  the  largest  spheres  to  have  the  slowest  uptake. 
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Figure  4.  Uptake  of  Anthracene  by  Paraffin  (partioning  basis). 


Figure  5.  Uptake  of  Anthracene  by  Paraffin  (adsorption  basis). 


The  data  for  anthracene  uptake  by  paraffin  on  an  adsorption  basis  (Figure  (5)) 
show  the  largest  paraffin  spheres  to  have  the  slowest  uptake.  This  is  as  expected  since 
anthracene  uptake  by  paraffin  is  not  on  an  adsorption  basis.  Similarly,  because  the  data 
represented  in  Figure  (3)  also  show  the  largest  nylon  66  spheres  to  have  the  slowest 
uptake,  anthracene  uptake  by  nylon  66  may  not  be  on  an  adsorption  basis. 

Since  the  sorption  mechanism  for  anthracene  by  nylon  66  could  not  be  defined 
with  certainty,  it  was  not  appropriate  to  fit  partitioning  parameters  to  the  data. 

Ceramic  Sphere  Rate  Study 

Solid  phase  concentrations  of  the  porous  ceramic  spheres  were  calculated  using 
Equation  (19): 

n 

-^10) 

- - - 7^ -  (19) 

"  5 


where: 

Vs  =  bulk  volume  of  ceramic  sphere. 

Experimental  data  was  plotted  as  (Cs  /Q)  versus  time  in  hours  as  seen  in  Figure  (6). 
Ceramic  and  Paraffin  Sphere  Rate  Study 

Solid  phase  concentrations  of  the  ceramic  and  paraffin  spheres  were  calculated 
using  (Equation  (20)): 


- 1 -Vkj)  )]  -Ci(n)Vi(„) 


Cs(n)  -  ■ 


H 


Vs  +  M, 


(20) 


Experimental  data  was  plotted  as  (Cs/Q)  versus  time  in  days  as  seen  in  Figure  (7). 
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Figure  6.  Uptake  of  Anthracene  by  Ceramic  Porous  Spheres. 


Figure  7.  Uptake  of  Anthracene  by  Ceramic  and  Paraffin  Spheres. 


The  expected  final  equilibrium  Cs(f)ICi(f)  for  the  porous  ceramic  sphere  rate  study 
was  the  porosity  of  T]  =  0.365.  Instead,  the  ratio  of  the  concentrations  approached  a  value 
greater  than  the  capacity  of  the  porosity  (Figure  (6)).  To  account  for  this  difference 
anthracene  must  not  only  be  diffusing  into  the  micropores  but  may  also  be  adsorbing  to 
the  silanol  groups  (SiOH)  of  the  kaolinite  mineral  structure.  This  behavior  is  well 
documented  and  has  been  measured  for  substituted  aromatic  compounds  on  kaolinite 
(Schwarzenbach  and  Westall,  1981).  If  the  process  for  this  unexpected  sorption  is 
adsorption,  additional  steps  to  prevent  the  adsorption  from  occurring  can  be  incorporated 
into  future  experimental  procedures.  The  silanol  groups  of  the  porous  ceramic  spheres 
can  be  deactivated  by  silanization  with  dimethylchlorosilane  (DMCS)  followed  by  a 
methanol  rinse.  This  reaction  replaces  the  silanol  hydrogen  with  dimethylmethoxysUane 
[SiOSi(CH3)20CH3]  which  minimizes  adsorption  by  solute  species.  Conversely,  one  may 
assume  that  the  additional  sorption  is  linear,  reversible,  and  fast  relative  to  the  slow, 
diffusive  transport  in  the  micropores  and  no  procedural  adjustments  are  necessary.  These 
assumptions  were  made  and  the  increased  sorption  was  included  in  the  MSS  model 
simulations. 

Blanks  and  Controls 

Instrumental  controls  were  analyzed  repeatedly  over  the  course  of  experimental 
analysis.  Controls  were  performed  before  experimental  samples  were  analyzed  and  again 
after  analysis  of  every  ten  experimental  samples.  The  experiments  that  were  analyzed  over 
the  course  of  many  weeks,  such  as  the  nylon  66  sphere  rate  study  and  the  mixed  ceramic 
and  paraffin  sphere  rate  study,  included  many  controls.  The  ceramic  sphere  rate  studies 
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were  conducted  over  the  course  of  hours  and  therefore  do  not  include  many  controls. 
With  the  exception  of  a  few  outliers,  instrument  controls  fall  within  a  95%  confidence 
interval. 

An  initial  condition  vial  was  sampled  repeatedly  over  the  course  of  the  nylon  66 
sphere  rate  study  to  determine  the  effects  of  loss  of  anthracene  due  to  volatilization, 
degradation,  and/or  sorption  on  vial  glass  and  septa.  After  69  days  a  loss  of  less  than  2% 
was  found.  Significant  volatilization  of  anthracene  is  unlikely  as  is  degradation  in  the 
methanol/water  solution.  Also,  the  anthracene  would  not  preferentially  sorb  to  the  glass 
from  the  liquid  50-/50  (v/v)  methanol/water  phase.  Loss  due  to  sorption  by  the  Teflon 
septa  is  possible,  yet  improbable.  The  amorphous  septa  would  possibly  swell  during 
exposure  to  the  liquid  phase  enhancing  sorptive  processes  as  in  lOMD  (Lion  et  al.,  1990). 
However,  during  the  experiment  the  vials  were  placed  upright  on  a  shaker  table  and  rarely 
inverted.  Solvent  contact  with  the  septa  was  minimal  and  not  of  length  to  induce  septa 
swell  and  subsequent  sorption.  The  loss  is  most  probably  due  to  operator  error  in  the 
experiment  or  eiror  in  instrument  analysis. 

Blanks  containing  the  solid  phase  of  study  and  anthracene  free  liquid  phase  were 
made  for  all  experiments.  In  all  cases  no  trace  of  anthracene  was  detected  in  the  blanks. 
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Model  Simulation  and  Validation: 


With  the  data  from  the  batch  sorption  experiments,  the  MSS  sorption  mass  transfer 
model  was  used  to  estimate  the  effective  diffusion  rate  coefficient  for  anthracene  in  the 
porous  ceramic  sphere  and  mixtures  of  porous  ceramic  and  paraffin.  The  numerical  code 
for  the  MSS  model  written  in  FORTRAN  (Appendix  C)  minimizes  the  weighted  sum  of 
square  for  error  residuals  (SSQ)  for  the  uptake  of  solute  by  using  the  liquid  phase 
concentration  (Equation  (19))  and  time  sampled.  Mass  transfer  parameters  were 
estimated  from  breakthrough  curve  data  by  performing  multiple  nonlinear  least  squares 
regressions  minimizing  SSQ  solving  for  the  parameters  that  best  fit  the  rate  data  (Meeter 
and  Wolfe,  1968).  The  model  converts  the  continuous  frequency  distribution  of  sorption 
sites  within  the  volume  of  the  sphere  (Equation  (1 1))  into  a  discrete  distribution.  Using 
synthetic  media,  the  model  will  optimally  converge  on  the  known  geometric  parameters 
(maximum  diffusion  path  length  and  shape  factor).  The  MSS  numerical  solution  was 
verified  using  an  analytical  solution  (Crank,  1975). 

The  MSS  model  was  used  to  estimate  the  effective  diffusion  coefficients  for  the 
homogeneous  experiments  involving  porous  ceramic  spheres.  The  model  estimated  the 
effective  diffusion  coefficients  using  the  known  maximum  diffusion  path  length,  dmu,  and 
shape  factor,  X,  for  the  porous  ceramic  spheres  and  an  average  Kp  found  in  independent 
experiments  with  the  porous  ceramic  spheres.  The  fitted  diffusion  coefficients  and  SSQ 
are  shown  in  Table  5. 
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Table  5.  Fit  Den  using  known  5,^  and  X. 


Solid  Phase 

Kp 

Deff(fit) 

(cmVsec) 

StdDev 

Smax 

(cm) 

■ 

SSQ 

0.55  cm  dia. 
(vial  1) 

0.425 

1.21E-07 

8.80E-09 

0.55 

2 

1.26E-03 

0.55  cm  dia. 
(vial  2) 

0.425 

2.25E-07 

3.45E-08 

0.55 

2 

1.40E-02 

Combined 
0.55  cm  dia. 

0.425 

1.65E-07 

1.20E-08 

0.55 

2 

1.67E-02 

0.75  cm  dia. 
(vial  1) 

0.425 

8.99E-07 

4.58E-08 

0.75 

2 

3.24E-03 

0.75  cm  dia. 
(vial  2) 

0.425 

9.29E-07 

8.71E-08 

0.75 

2 

1.05E-02 

Combined 
0.75  cm  dia. 

0.425 

9.19E-07 

4.82E-08 

0.75 

2 

1.35E-02 

Considering  the  model  results  (Table  (5))  and  contrary  to  that  anticipated,  it 
appears  that  diffusion  in  the  smaller  spheres  (0.55  cm  dia.)  is  occurring  at  a  faster  rate  than 
for  the  larger  spheres  (0.75  cm  dia).  This  is  possibly  an  artifact  of  the  incomplete 
experiment  small  sphere  experiment  as  unfortunately  the  small  sphere  data  did  not  reach 
equilibrium  before  the  HPLC  suffered  mechanical  failure.  However,  using  the  Kp  from  the 
large  spheres  to  fit  the  should  have  accounted  for  the  lack  of  a  complete  data  set  and 
resulted  in  a  similar  for  the  two  experiments.  It  is  possible  that  the  Kp  for  the  small 
spheres  is  not  the  same  as  the  Kp  for  the  large  spheres  although  unlikely.  A  more  probable 
reason  for  the  difference  in  Dg^is  in  the  experimental  conditions.  Prelinoinary  experimental 
mns  with  the  solute  and  sorbent  using  shaker  tables  to  establish  well-mixed  conditions 
resulted  in  a  cloudy  liquid  phase  that  clogged  HPLC  tubing  and  filters.  Subsequent  runs 
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were  not  agitated,  merely  swirled.  It  is  possible  the  difference  in  Deff  between  the  large 
and  small  spheres  is  due  to  insufficient  mixing. 

The  results  for  Deff  in  the  large  spheres  (0.75  cm  dia.)  were  compared  with  the 
fitted  diffusion  coefficients  of  previous  research  using  the  same  size  porous  ceramic  media, 
and  solutes,  modeled,  however,  with  a  different  analytical  model  (Rao  et  al., 
1980).  Before  a  direct  comparison  can  be  made,  the  Deff  for  the  anthracene  (a  sorbing 
solute)  was  adjusted  to  account  for  retardation  due  to  sorption  when  compared  to  the 
nonsorbing  solutes  used  by  Rao  et  al.  (1980).  This  adjustment  was  made  with  Equation 
(21)  where  the  effective  diffusion  for  the  sorbing  solute  is  multiplied  by  the  retardation 
factor,  R,  for  the  sorbing  solute  to  give  an  equivalent  effective  diffusion  as  though  the 
solute  were  nonsorbing. 

^efflnonsorbing)  ~  (sorbing)  (21) 

The  diffusion  coefficient  found  in  this  effort  is  an  order  of  magnitude  smaller 
(Table  6.)  than  that  found  by  Rao  et  al.  (1980).  This  may  indicate  that  diffusive  transport 
in  this  experiment  was  somewhat  slower  than  for  the  work  by  Rao  et  al.  (1980b),  possibly 
due  to  the  differences  in  solutes  for  the  experiments. 


Table  6.  Effective  Diffusion  Coefficients  for  Nonsorbing  Solutes  into  0.75  cm  dia.  Porous 
Ceramic  Spheres. 


Diffusion  Coefficient 
cm^/hr 

Anthracene 

36q. 

'H2O 

Deff 

0.0039 

0.012 

0.012 

Do 

0.0316 

0.0396 

0.0396 
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The  tortuosity,  t,  is  represented  by  a  ratio  of  the  effective  diffusion  to  the  diffusion 


coefficient  for  a  given  media  (Equation  (22)). 


T  = 


(22) 


Investigating  the  ratio  of  the  effective  diffusion  to  the  diffusion  coefficient  for  each 
experiment  respectively  results  in  a  x  for  this  experiment  equal  to  0.122  and  for  the  work 
by  Rao  et  al.  (1980)  equal  to  0.303.  It  would  be  expected  that  for  the  systems  using  the 
same  sorbent,  x  would  be  equivalent.  Although  relatively  close,  the  values  are  not 
equivalent.  One  possible  reason  for  the  variation  may  be  due  to  steric  hindrance  due  to 
the  comparatively  large  size  of  the  anthracene  molecules.  However,  because  the  size  of 
the  “micropores’  of  the  ceramic  spheres  are  presumed  large  relative  to  the  anthracene 
molecules,  this  would  be  unlikely.  Another  possible  and  perhaps  more  probable  reason 
could  be  due  to  experimental  technique.  In  both  this  experiment  and  the  Rao  et  al.  (1980) 
experiment  the  vials  were  swirled  intermittently.  Perhaps  the  swirling  techniques 
employed  in  this  experiment  were  not  sufficient  to  establish  a  well-mixed  liquid  phase,  thus 
impacting  diffusive  transport. 

For  real  soils,  the  media  shape  properties  are  unknown.  Real  soil  is  a 
heterogeneous  mix  of  particles  with  differing  partitioning  coefficients,  Kp,  effective 
diffusion,  Deff,  and  distribution  of  sorption  sites,/(5).  The  model  must  be  able  to 
overcome  these  unknowns  by  using  sorption  characteristics  that  is  readily  available  or 
assumed,  and  fitting  as  few  remaining  soil  parameters  as  possible.  Using  an  average  Kp  for 
the  soil  and  the  free  liquid  diffusivity  for  the  solute  of  interest,  the  model  can  be  used  to  fit 
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parameters  of  the  average  soil  particle.  This  average  or  virtual  particle  can  be  defined  by 
the  shape  properties  of  a  virtual  maximum  path  length,  5*ma,  and  some  shape  factor,  A,. 
Converting  the  maximum  path  length  is  necessary  to  adjust  for  model  use  of  the  free  liquid 
diffusivity  to  fit  the  soil  parameters.  The  real  maximum  path  length,  8max,  was  converted 
to  a  virtual  maximum  path  length,  5*^,  using  Equation  (23): 

^  max  ~  ^  max  (23) 

Do  is  the  free  liquid  diffusivity  of  anthracene  in  50/50  (v/v)  methanol/water  which 
was  used  in  place  of  the  unknown  effective  diffusion  coefficient  to  fit  the  shape 
parameters.  Table  7  shows  the  MSS  model  results  using  the  free  liquid  diffusivity  to  fit 
S' max  and  X  for  the  heterogeneous  system  of  mixed  paraffin  and  porous  ceramic  spheres. 


Table  7.  Two-parameter  fit,  6  and  X,  using  free-liquid  diffusion  coefficient.  Dp. 


Solid  Phase 

Kp 

Do 

(cm^/sec) 

C* 

O  max 

(fit) 

(cm) 

■ 

X  (fit) 

SSQ 

PMB  1 

0.888 

8.37E-06 

159.4 

25.2 

60.1 

14.6 

1.45E-01 

PMB2 

0.888 

8.37E-06 

87.4 

16.8 

31.6 

7.9 

1.39E-02 

Combined  PMB 

0.888 

8.37E-06 

122.4 

16.5 

46.5 

8.75 

4.07E-02 

Because  the  actual  media  used  is  of  known  geometric  properties,  the  true 
distribution  of  sorption  sites  is  also  known.  The  graphical  representation  of  the  tme 
distribution  of  sorption  sites  was  developed  using  a  composite  of  the  sorbent  geometry  for 
each  media  in  the  system;  0.55  cm  diameter  porous  ceramic  spheres,  0.75  cm  diameter 
porous  spheres  and  0.499  cm  diameter  paraffin  spheres.  The  mass  fraction  contribution  of 
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each  media  type  were  summed  to  create  a  resultant  combined,  true  distribution  (Figure 

(8)). 

Using  a  distribution  of  sorption  capacity  along  the  path  length  as  defined  by 
Equation  (11),  the  fitted  values  of  5*max  and  A,  were  plotted  and  compared  to  the  tme 
average  distribution  of  sorption  sites  (Figure  (9)).  The  model  fit  is  not  coincident  with 
that  of  the  true  distribution  for  the  media  except  at  long  path  lengths.  One  possible  reason 
for  this  is  in  the  computation  of  the  true  distribution  of  sorption  capacity  sites  along  the 
path  length.  The  assumptions  entered  into  in  this  investigation  and  deemed  necessary  for 
analysis  centered  on  the  use  of  a  synthetic  media.  By  using  synthetic  media  it  was 
assumed  that  soil  sorption  characteristics  and  geometric  characteristics  were  known.  It  is 
possible  that  the  experiment  as  conducted  was  insufficient  to  establish  known  soil  sorption 
characteristics.  For  example,  the  experiment  was  unable  to  be  run  until  equilibrium 
conditions  were  reached.  Also,  as  was  previously  suggested,  it  is  possible  that  conditions 
were  not  sufficient  to  establish  a  well-mixed  liquid  phase.  Thus,  because  the  “true” 
distribution  may  have  been  developed  with  inaccurate  data,  it  cannot  be  said  that  the  soil 
sorption  characteristics  are  known.  If  this  is  the  case,  the  model  has  cannot  be  evaluated 
as  there  is  no  known  true  distribution  of  sorption  sites.  Future  experiments  to  assess  the 
model  validity  must  ensure  experimental  conditions  are  such  that  soil  sorption 
characteristics  are  accurate  so  that  they  may  be  used  to  develop  a  true  distribution  of 
sorption  sites. 
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Figure  8.  True  “Virtual”  Geometry  created  from  Composite  of  Paraffin  and  Porous  Ceramic  Spheres. 


Distributions  of  Sorption  Capacity  for  Mixed  Media  Experiment. 


V.  Summary  and  Conclusion 


The  sorption  rates  of  anthracene  on  synthetic  media  consisting  of  nylon  66,  porous 
ceramic,  paraffin  and  mixtures  of  porous  ceramic  and  paraffin  spheres  were  measured 
using  HPLC.  The  results  were  interpreted  by  use  of  a  multisite  diffusion  mass  transfer 
model. 

The  time  required  to  establish  equilibrium  conditions  for  the  nylon  66  was  in 
excess  of  the  allotted  time  for  this  research  effort.  The  low  rates  of  uptake  suggest  the 
nylon  66  to  be  a  highly  condensed  polymer  where  the  time  required  to  reach  equilibrium  is 
much  greater  than  for  less  condensed  polymers  such  as  the  paraffin. 

The  rate  behavior  in  the  porous  ceramic  spheres  was  slower  than  published  results 
using  the  same  media  although  this  may  be  attributed  to  experimental  techniques.  The  rate 
parameters  infenred  from  the  data  were  not  entirely  consistent  among  sizes  of  spheres, 
probably  attributed  to  experimental  deficiencies  where  equilibrium  conditions  were  not 
sampled  due  to  mechanical  failure  of  the  analytical  instrumentation. 

These  results  demonstrate  the  significance  of  adequate  experimental  conditions  to 
establish  accurate  depiction  of  the  distribution  of  sorption  sites  in  the  synthetic  media  to 
allow  model  comparison.  The  experimental  results  reported  here  have  important 
implications  for  the  proper  design  of  corresponding  laboratory  experiments. 
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Appendix  A:  Glossary 


Roman 


b  =  radius  of  sphere,  L 

C  =  concentration  of  individual  components,  ML'^ 

Cim  =  concentration  (path  length  5  and  time  t)  in  the  immobile  liquid  phase,  ML^ 

Cim  =  average  concentration  in  the  immobile  liquid  phase,  ML'^ 

Cl  =  concentration  of  solute  in  solution,  ML'^ 

Ci(i)  =  initial  liquid  phase  concentration,  ML'^ 

Cm  =  concentration  in  the  mobile  liquid  phase,  ML^ 

Com  =  concentration  ofsorbate  associated  with  the  SOM,  ML'^ 

C,(j)  = final  liquid  phase  concentration,  M  L'^ 

^  =  diffusion  coeffecient,  t'^ 

Da  -  apparent  diffusion  coefficient,  ll  t^ 

Dm  =  mobile  region  dispersion  coefficient,  L?  t^ 

Do  =  free-liquid  diffusion  coefficient,  l}  t^ 

f= fraction  of  sorption  sites  that  are  in  contact  with  the  mobile  water  phase 

fee  =  fraction  of  organic  carbon 

fom  =  weight  fraction  of  solid  which  is  SOM 

F  =  fraction  of  sorption  sites 

J  =  mass  flux,  M  L'^  f‘ 

ki  =  first-order  mass  transfer  rate  coefficient,  L  t^ 

Kad  =  equilibrium  adsorption  coefficient,  M'^ 

Kd  =  linear  sorption  distribution  coefficient,  lI  M'^ 

Koc  =  organic  carbon  distribution  coefficient,  M'‘ 

Kp  =  equilibrium  partition  coefficient,  M'^ 

Kow  =  octanol-water  partition  coefficient,  if  M'‘ 
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Ml  =  mass  of  the  liquid,  M 
Ms  =  mass  of  the  solid,  M 
r  =  radial  coordinate  of  the  spherical  particle,  L 
t  =  time,  t 

R  =  retardation  factor 
S  =  sorbed  concentration,  M  M'^ 

V  =  average  pore  water  velocity,  L 
Vs  =  bulk  volume  of  ceramic  sphere 
Vi  =  volume  of  the  liquid, 

X  =  distance,  L 

Greek 

a '  =  first-order  mass  tranter  coefficient,  L 

5  =  diffusion  path  length,  L 

S/nax  *  maximum  diffusion  path  length,  L 

A,  =  shape  factor 

p*  =  bulk  soil  density,  ML'^ 

p:  =  density  of  the  liquid,  M  L'^ 

0m  =  mobile  phase  porosity,  if  L'^ 

0im  =  immobile  phase  porosity,  L'^ 

X  =  tortuosity 
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Appendix  B:  Experiment  Data 


Table  1.  Nylon-66  Isotherm  Experiment 


Liquid  Phase  Concentration  of 
Anthracene 

Contents  of  each  Vial 
~  0.07  g  crushed  Nylon  66 

0.81  mg/L 

BN-2A 

BN-2B 

BN-5A 

BN-5B 

0.61  mg/L 

BN-2A 

BN-2B 

BN-5A 

BN-5B 

0.41  mg/L 

BN-2A 

BN-2B 

BN-5A 

BN-5B 

0.20  mg/L 

BN-2A 

BN-2B 

BN-5A 

BN-5B 

0.0  mg/L 

BN-2 

BN-5 

Table  2.  Nylon-66  Rate  Experiment  (BN1.5  -BN5) 


Size  of 

Nylon  66  Sphere 

(cm) 

Quantity  of 
Spheres  per  Vial 

Initial 

Concentration 
of  Liquid  Phase 
(mg/L) 

Quantity  of 
Vials 

Mass  of 
Solid  Phase 

(g) 

0.238 

71 

0.81 

3 

0.572 

0.318 

30 

0.81 

3 

0.578 

0.476 

9 

0.81 

3 

0.560 

0.635 

4 

0.81 

3 

0.612 

0.794 

2 

0.81 

3 

0.596 

mixed  sizes 

3  (0.238),  2  (0.318), 

1  (0.476),  1  (0.635), 

1  (0.794) 

0.81 

3 

0.576 

mixed  blank 

3  (0.238),  2  (0.318), 

1  (0.476),  1  (0.635), 

1  (0.794) 

0.81 

2 

0.576 

Table  3.  Porous  Ceramic  Sphere  Rate  Experiment  (SB  and  LB) 


Size  of 

Porous  Ceramic 
Sphere 
(cm) 

Quantity  of 
Spheres  per  Vial 

Initial 

Concentration  of 
Liquid  Phase 
(mg/L) 

Quantity 
of  Vials  or 
Bottles 

Volume  of 
Liquid  Phase 

(mL) 

0.55 

8 

0.80 

2 

19.0 

0.75 

6 

0.80 

2 

35.0 
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Table  4.  Mixed  Porous  Ceramic  and  Paraffin  Sphere  Rate  Experiment  (PMB) 


Size  of  Sphere 

(cm) 

Quantity  of 
Spheres  per  Vial 

Initial 

Concentration  of 
Liquid  Phase 
(mg/L) 

Quantity 
of  Vials 

Volume  of 
Liquid  Phase 

(mL) 

0.55,  0.75, 

6  (0.55),  3  (0.75), 

ceramic  and 
0.249  paraffin 

and  7  (0.249) 

0.80 

2 

31.0 
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Appendix  C:  MSS  Model  Code 


MSSHT 

C 

C  Program  Testfit 
C 

C  Adapted  for  CSTR  experiencing: 

C  Multisite  (series)  sorption  on  soil  (spherical  distribution) 

C  Flow  and  inlet  concentration  perturbations 

C  Bicontinuum  sorption  on  reactor  walls 

C 

C  NUMERICAL  SOLUTION 
C 

C  BY:E.  Heyse 

C  University  of  Florida,  23  AUG  93 

C 

C 

C  dim  scrat  =  5*np  +np*np  +(2+np)*nob 
C  signs  =  1  means  parms  cannot  change  sign 
C  P(IPOINT(I))  =  VAL(I)  ==  PNAME(I) 

C 

C  IMPLICIT  DOUBLE  PRECISION  (A-H,aZ) 

DOUBLE  PRECISION  FOC,XM,V,PC0,TBCC(30),XQ(30),XCIN(30) 

DOUBLE  PRECISION  VR(30),VA(30),CA(30)J^AST 

DOUBLE  PRECISION  XMR,XKPR,XKR,XFR,X(100),XKP,DEFF,DPATH,XSS0,XSR0 
DIMENSION  IX(100),Y(100),P(7),DIFF(7),SIGNS(7),SCRAT(1999) 

DIMENSION  IFC(30) 

CHARACTER*8  PNAME(7) 

CHARACTER*  15  INFIL,OUTl,OUT2 
DIMENSION  VAL(7),IPOINT(7),IVARY(7) 

COMMON /BMOD/FOC,XM,V,PCO,NEFD,TBCC,XQ,XCIN,VR,VA,CA 

COMMON  /BMODA/XMR,XKPR,XKR,XFR,NFRAC,IFLAG 

COMMON  /VALUE/  VAL,IPOINT,IVARY,PNAME,NnX 

COMMON  /BDAT/  IX, X 

CHARACTER*80  TIT,LIN1,LIN2,LIN3,LIN4 

NFIX=6 

PNAME(1)='  Kd ' 

PNAME(2)='Deff' 

PNAME(3)='dpath' 

PNAME(4)='  Lambda ' 

PNAME(5)='  SsO ' 

PNAME(6)='  Fast ' 

WRITE(*,*)  'NAME  OF  INPUT  FILE?' 

READ(5,1111)INFIL 
nil  FORMAT(A15) 

WR1TE(6,*)'NAME  OF  SHORT  OUTPUT  HLE?' 

READ(5,llll)OUTl 
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WRITE(6,*)'NAME  OF  LONG  OUTPUT  FILE?' 
READ(5,llll)OU'r2 

OPEN(UNIT=3,FILE=INHL,STATUS=’OLD’) 
OPEN(UN1T=10,FILE=OUT1  ,STATUS='NEW’) 
OPEN(UNIT=4,nLE=OUT2,STATUS=’NEW’) 

READ(3;(A)')  "HT 
C 

C  Read  parameters  to  be  fitted 
C 

C  RKD  =  partition  coefficient,  KP, 

C  DEFF  =  Coefficient  of  Diffusion,  D_eff  cm'^2/sec 
C  DPA'TH  =  Total  diff.  path  length,  d,  cm 
C  RLAM  =  shape  factor  Lambda 

C  XSSO  =  Initial  dmls  Concentration  in  slow  soil  sites,  Sso,  ug/g 
C  FAST  =  Fraction  of  equilibrium  sorption  sites 
C 

READ(3,’(A)')  LINl 

READ(3,*)  RKD,DEFF,DPA'rH,RLAM,XSSO  J^AST 
C 

C  Read  Constants 
C 

C  FOC  =  Mass  fraction  organic  carbon  on  soil,  foe 
C  XM  =  Mass  of  soil  in  reactor,  g 
C  V  =  Volume  of  solvent  in  reactor  at  time  zero,  mL 
C  PCO  =  initial  concentration  in  liquid,  mg/1 
C  NEFD  =  Number  of  boundary  condition  changes 
C  EFC  =  boundary  condition  number 
C  TBCC  =  time  for  boundary  condition  to  end 
C  XQ  =  Flow  rate,  mLAnin 
C  XCIN  =  Dimensionless  inlet  concentration 
C  XMR  =  Mass  of  reactor,  g 
C  XKPR  =  Linear  partition  coef.  on  reactor,  mL/g 
C  XKR  =  1st  order  mass  transfer  rate  onto  reactor,  min'^-1 
C  XFR  =  Fraction  instaneous  sorbing  reactor  sites 
C  VS(J)  =  volume  of  liquid  in  rx  for  BC  period  j,  ml 
C  VR(J)  =  volume  of  liquid  removed  at  end  of  BC  period  J,  ml. 

C  VA(j)  =  volume  of  liquid  added  at  end  of  BC  period  J,  ml. 

C  CA(J)  =  concentration  in  liquid  added  at  end  of  BC  period  J ang/l. 
C  XIC  =  sorbed  concentration  in  each  slow  compartment,  ug/g 
C  CL  =  cument  liquid  concentration,  mgA 
C 

READ(3,'(A)')  LIN2 
READ(3,*)  (IVARY(I),I=1,NHX) 

READ(3,'(A)')  LIN3 
READ(3,*)  XM,V,PC0 
C  READ(3,*)  XMR,XKPR,XKR,XFR 
READ(3,*)  NEFD 
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D0  6I=1,NEFD 

6  READ(3,*)  IFC(I),TBCC(I),VR(I),VA(I).CA(I) 

C  Read  numerical  parameter 
C  IFLAG=0  for  sorbed/liquid  concentration 
C  IFLAG=1  for  liquid  concentration 
READ(3,*)NFRAC,IFLAG 
C 

C  Read  Data  to  be  fitted 
C 

C  X=  time,  seconds 
C  Y  =  observation  at  time  X,  see  IFLAG. 

C 

READ(3,'(A)')LIN4 

NOB=0 

DO  3  1=1,100 

READ(3,*JEND=4)  X(I),Ya) 

3  NOB=NOB+l 

4  WRITE(4,1)  TIT 

1  FORMATC/'  ’A) 

WRITE(4,1)LIN1 

WRITE(4,2)  RKD,DEFF,DPATH,RLAM,XSSO,FAST 

2  FORMATC  ',7E12.6) 

WRITE(4,1)LIN2 
WRnE(4,5)  (IVARY(D,I=1JSHX) 

5  FORMATC  ’,7110) 

WRITE(4,1)  LIN3 
WRITE(4,777)  XM,V,PC0,NEFD 

777  FORMATC  '.(3E12.6,/),I10) 

DO  778I=1,NEFD 

WRnE(4,779)  IFC(I),TBCCa),VR(I),VA(I),CA(I) 

778  CONTINUE 

779  FORMATC  ',I10,4E12.6) 

VAL(1)=RKD 

VAL(2)=DEFF 

VAL(3)=DPATH 

VAL(4)=RLAM 

VAL(5)=XSS0 

VAL(6)=FAST 

NP=0 

JHX=NnX 
DO  20 1=1, NHX 
IF  aVARY(I).EQ.l)  THEN 
NP=NP+1 
IPOINT(I)=NP 
ELSE 

IPOINT(I)=JFIX 

JFIX=JFIX-1 
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ENDEF 

20  CONTINUE 
DO  22  I=1,NFIX 
P(rPOINT(I))=VAL(I) 

22  CONTINUE 
NPROB=I 
MIT=0 
EPS1=0.0 
EPS2=0.0 
FLAM=0.0 

IF(EPSl.LT.l.E-9)  EPS1=0.0 
IF(EPS2.LT.I.E-9)  EPS2=I.E-3 
IF(MIT.EQ.O)  MIT=I5 
IF(FLAM.LT.l.E-9)  FLAM=0.01 
FNU=10. 

DO  9  1=1,6 
SIGNS(I)=1.0 
9  DIFF(I)=0.01 

WRnE(4,7)  NPROB 

7  FORMATC  PROBLEM  NUMBERM3,T26,’X’,T38,'Y’/) 

WRnE(4,8)  (X(I).Y(I),I=1NI0B) 

8  FORMATC  ’,T19,2E12.5) 

WRITE(*,12)  (PNAMEa),I=l  WX) 

WRrrE(10,12)  (PNAME(D,I=1.NFIX) 

12  FORMAT(yT6,'SSQ',T11.7A12) 

WRITE(*,15)  (P(IP0INT(D),I=12^FIX) 

WRITE(10,15)  (P(IPOINTa)),I=l,NnX) 

15  F0RMAT(T11,7E12.5) 

CALL  UWHAUS(NPROB,NOB,YJ^PTDIFF,SIGNSJEPSl,EPS2>nTELAM, 
*FNU,SCRAT) 

WRITE(*,14)  TIT 
WRITE(10,14)  TIT 
14  FORMATC/"  END:  ',A) 

999  STOP 
END 

SUBROUTINE  UWHAUS(NPROB  J^OB,Y.NP,TH, DIFF, SIGNS, EPS  1  ,EPS2, 

1  MIT,  FLAM,  FNU,  SCRAT) 

DIMENSION  SCRAT(1),Y(1),TH(1),DIFF(1),SIGNS(1) 

IA=1 

IB=IA+NP 

IC=IB+NP 

ID=IC+NP 

IE=ID+NP 

IF=IE4-NP 

IG=IF+NOB 

IH=IG+NOB 
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II  =  IH  +  NP*NOB 
IJ  =  IH 

CALL  HAUS59(NPROB,NOB,Y.NP,TH,DIFF, SIGNS, EPS1,EPS2MIT 

1  JFLAM,FNU,SCRAT(IA),  SCRAT(IB),  SCRAT(IC),  SCRAT(ID), 

2  SCRAT(IE),  SCRAT(IF),  SCRAT(IG).  SCRAT(IH),  SCRAT(II), 

3  SCRAT(IJ)  ) 

RETURN 

END 

SUBROUTINE  HAUS59(NPRBO,  NBO,  Y,NQ,TH,DIFZ,SIGNS,EP1S,EP2S, 
1MIT,FLAM,FNU,  Q,P,EJ>HLTBE,R,A,D,DELZ) 

DIMENSION  TH(1),  DIFZ(l).  SIGNS(l),  Y(l),  Q(l),  P(l),  E(l), 

1  PHI(l),  TB(1),  F(l),  R(l).  A(l),  D(l),  DELZ(1) 

CHARACTER*8  PNAME(7) 

DIMENSION  VAL(7),IPOINT(7),IVARY(7) 

COMMON  A^ALUE/  VAL,IPOINT,IVARY,PNAME,NHX 
DATAMAXCNT/6/ 

ACOS(X)  =  ATAN(SQRT(L0/X**2  -  1.0)) 

NP  =  NQ 

NPROB  =  NPRBO 
NOB  =  NBO 
EPS1  =  EP1S 
EPS2  =  EP2S 
NPSQ  =  NP*NP 

NSCRAC  =  5*NP+NPSQ  +2*NOB+NP*NOB 
WRITE(4,  1000)  NPROB,  NOB,  NP,  NSCRAC 
WRITE(4,  1001) 

CALL  GASS60(1,  NP,  TH,  TEMP,  TMEP) 

WRITE(4,  1002) 

CALL  GASS60(1,  NP,  DIFZ,  TEMP,  TEMP) 

IF(MIN0(NP- 1 ,50-NP,NOB-NP,M[T- 1 ,999-MIT))99, 15, 15 

15  IF(FNU-1.0)99,99, 16 

16  CONTINUE 
DO  19  I=1,NP 
TEMP  =  ABS(DIFZ(I)) 

IF(AMIN1(  1.0-TEMP,  ABS(TH(I))))99, 99, 19 
19  CONTINUE 
GA  =  FLAM 
NIT=  1 

ASSIGN  131  TO  LAOS 
ASSIGN  225  TO  IRAN 
ASSIGN  265  TO  JORDAN 
IF(EPS1)  5,  10,  10 
5  EPSl  =  0 
10  IF(EPS2)  40,  40,  30 
40  IF(EPS1)60,  60,  50 
60  ASSIGN  270  TO  IRAN 
GO  TO  70 
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50  ASSIGN  265  TO  IRAN 
GO  TO  70 

30  IF(EPSl)  80,  80,  70 
80  ASSIGN  270  TO  JORDAN 
70  SSQ  =  0 
D0  71I=l,NnX 

IFaVARYa).EQ.l)  VAL(I)=TH(IPOINT(I)) 

71  CONTINUE 

CALL  MODEL(NPROB,  TO,  F,  NOB,  NP) 

DO  901=  l,NOB 
R(I)  =  Y(D-Fa) 

90  SSQ=SSQ+Ra)*R(I) 

WRITE(4,  1003)  SSQ 

WRITE(*,3003)  SSQ,(TO(IPOINT(I)),I=l,NFIX) 
WRITE(10,3003)  SSQ.(TH(IPOINT(I)),I=l,NFIX) 

3003  FORMATC  ',lPE9.2,0P7El'2.5) 

WRITE(4,  1011) 

WRITE(4,  2006)  (F(I),I=1  JJOB) 

C  BEGIN  ITERATION 

100GA  =  GA/FNU 
INTCNT  =  0 
WRITE(4,  1004)  NIT 
101  JS  =  1  -  NOB 
DO  130J=1,NP 
TEMP  =  TH(J) 

P(J)=DIFZ(J)*TO(J) 

TO(J)=TH(J)+P(J) 

Q(J)=0 

JS  =  JS  -H  NOB 
DO  102I=1,NFIX 

IF(IVARY(I).EQ.l)  VAL(I)=TO(IPOINT(I)) 

102  CONTINUE 

CALL  MODEL(NPROB,  TO,  DELZ(JS),  NOB,  NP) 

IJ  =  JS-1 

DO  120 1=1,  NOB 
U  =  IJ  +  1 

DELZ(U)  =  DELZ(U)  -  F® 

120  Q(J)  =  Q(J)  +  DELZ(U)  *  R® 

Q(J)=Q(J)/P(J) 

C  Q=XT*R  (STEEPEST  DESCENT) 

130TO(J)  =  TEMP 
GO  TOLAOS,(131,414) 

131  DO150I=l,NP 
DO  151  J=1,I 
SUM  =  0 
KJ  =  NOB*(J-l) 

KI  =  NOB*(I-l) 
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DO  160K=  l,NOB 
KI  =  KI  +  1 
KJ  =  KJ  +  1 

160  SUM  =  SUM  +  DELZ(KI)  *  DELZ(KJ) 

TEMP=  SUM/(P(I)*P(J)) 

JI  =  J  +  NP*(I-1) 

D(JI)  =  TEMP 
IJ  =  I  +  NP*(J-1) 

151  Dai)  =  TEMP 
150  E®  =  SQRT(D(JI)) 

666  CONTINUE 
DO  153  1=  1,NP 
IJ  =  I-NP 
DO  153J=1,I 
IJ  =  IJ  +  NP 

A(IJ)  =  D(U)/(E(D*E(J)) 
n  =  J  +  NP*(I-l) 

153  A(JI)  =  A(IJ) 

C  A=  SCALED  MOMENT  MATRIX 

II  =  -NP 
DO  155  1=1, NP 
P®=Qffi/E® 

PHI(I)=P(I) 

II  =  NP+  1  +  n 
155  A(II)  =  A(II)  +  GA 
1=1 

CALL  MATTNCA,  NP,  P,  I,  DET) 

C  P/E  =  CORRECTION  VECTOR 

STEP=1.0 
SUM1=0. 

SUM2=0. 

SUM3=0. 

DO  231 1=1, NP 
SUM1=P(I)*PHI(I)+SUM1 
SUM2=P(I)*P®+SUM2 
SUM3=  Pffla)  *  PHI(I)  +  SUM3 
231  PHI(I)  =  Pa) 

TEMP  =  SUM1/SQRT(SUM2*SUM3) 

TEMP  =  AMINUTEMP,  1.0) 

TEMP  =  57.295*ACOS(TEMP) 

WRITE(4,  1041)  DET,  TEMP 
170  DO220I=  1,NP 
P(I)  =  PHI(I)  *STEP/E(I) 

TB(I)  =  TH(I)  +  P(I) 

220  CONTINUE 
WRITE(4,  7000) 

7000  FORMATC  TEST  POINT  PARAMETER  VALUES’) 
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WRITE(4,  2006)  (TB(I),  I  =  1,  NP) 

DO  221  1=  1,NP 
IF(SIGNS(I))  221,  221,  222 

222  IF(SIGN(1.0,TH(I))*SIGN(1.0,TB(I)))  663, 221,  221 
221  CONTINUE 

SUMB=0 
DO  223  1=1, NHX 

IF(IVARY(I).EQ.l)  VAL(I)=TB(IPOINT(I)) 

223  CONTINUE 

CALL  MODEL(NPROB,  TB,  F,  NOB,  NP) 

DO  230 1=1, NOB 
Ra)=Ya)-F(i) 

230  SUMB=SUMB+R(I)*Ra) 

WRITE(4,  1043)  SUMB 
WRITE(*,3002)  SUMB,(VAL(I),I=1OTX) 
WRITE(10,3002)  SUMB,(VALa),I=l,NFIX) 

3002  FORMATC  ',1PE9.2,0P7E12.5) 

WRITE(4,  1011) 

WRITE(4,  2006)  (F(I),I=l,NOB) 

IF(SUMB  -  (L0+-EPS1)*SSQ)  662, 662,  663 

663  IF(AMIN1(TEMP-30.0,GA))665,665,664 
665  STEP=STEP/2.0 

INTCNT  =  INTCNT  +  1 
IF(INTCNT-MAXCNT)  170,2700,2700 

664  GA=GA*FNU 
INTCNT  =  INTCNT  +  1 
IF(INTCNT-MAXCNT)  666,2700,2700 

662  WRITE(4, 1007) 

DO  669I=1,NP 
669  THa)=TB(I) 

CALL  GASS60(1,  NP,  TH,  TEMP,  TEMP) 

WRITE(4,  1040)  GA,  SUMB 
GO  TO  IRAN,(225,270,265) 

225  DO240I=  1,NP 

IF(ABS(P(I))/(l.E-20+ABS(TH(I)))-EPS2)  240,  240,  241 
241  GO  TO  JORDAN,(265,270) 

240  CONTINUE 
WRITE(4,  1009)  EPS2 
GOTO  280 

265  IF(ABS(SUMB  -  SSQ)  -  EPS1*SSQ)  266,  266, 270 

266  WRITE(4,  1010)  EPSl 
GO  TO  280 

270  SSQ=SUMB 
NIT=NIT+1 

IF(NIT  -  MIT)  100,  100,  280 
2700  WRITE(4,  2710) 
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2710  FORMAT(//l  15H  ****  THE  SUM  OF  SQUARES  CANNOT  BE  REDUCED  TO  THE 
SUM 

I  OF  SQUARES  AT  THE  END  OF  THE  LAST  ITERATION  -  ITERATING  STOPS  /) 

C  END  ITERATION 

280  WRITE(4,  1011) 

WRITE(4,  2001)  (F(I),  I  =  1,  NOB) 

WRITE(4,  1012) 

WRITE(4,  2001)  (R(I),  I  =  1,  NOB) 

SSQ=SUMB 

IDF=NOB-N? 

WRITE(4,  1015) 

1=0 

CALL  MATIN(D,  NP,  P,  I,  DET) 

DO  7692 1=1  NJP 

II  =  I  +  NP*(I-1) 

7692  E(I)  =  SQRT(Dai)) 

DO  340I=1,NP 
JI  =  I  +  NP*(H)  - 1 
IJ  =  I  +  NP*(I-2) 

DO  340  J  =  I,  NP 
JI  =  n  +  1 

A(n)  =  D(JD/(E(I)*E(J)) 

IJ  =  I J  +  NP 

340  AaJ)  =  A(JI) 

CALL  GASS60(3,  NP,  TEMP,  TEMP,  A) 

WRITE(4,  1016) 

CALL  GASS60(1,  NP,  E,  TEMP,  TEMP) 

IF(IDF)  341,  410,  341 

341  SDEV  =  SSQ/IDF 
WRITE(4,  1014)  SDEV,  IDF 
SDEV  =  SQRT(SDEV) 

WRITE(4,3007)  (PNAME(I),I=1,NHX) 

WRITE(*,3007)  (PNAMEa),I=l,NHX) 

WRITE(10,3007)  (PNAME(I),I=1,NFIX) 

3007  FORMAT(yT6,'SSQ',Tll,7A12) 

DO  3008  1=1, NHX 

IFaVARY(I).EQ.l)  VALa)=TH(IPOINTa)) 

3008  CONTINUE 

WRITE(*,3004)  SSQ,(VALa),I=l,NHX) 

WRITE(10,3004)  SSQ,(VAL(I),I=1,NFIX) 

WRITE(4,3004)  SSQ,(VAL(I),I=1,NHX) 

3004  FORMAT(lPE10.2,0P7E12.5) 

DO  3005  1=1, NP 

3005  P(I)=E(I)*SDEV 
DO  3009  I=1,NFIX 
VAL(I)=0.0 

IF(IVARY(I).EQ.l)  VAL(I)=P(IPOINTa)) 
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3009  CONTINUE 

WRITE(*.3006)  (VAL(I),I=1,NFIX) 

WRIIEC  10,3006)  (VAL(I), 1=1  NOFIX) 

WRITE(4,3006)  (VALa),I=  UNFIX) 

3006  FORMATC  Std.  Dev.',7E12.5) 

DO  391 1=1, NP 
P(I)=TH(I)+2.0*E(I)*SDEV 
391  TB(I)=TH(I)-2.0*E(I)*SDEV 
WRITE(4,  1039) 

CALL  GASS60(2,  NP,  TB,  P,  TEMP) 

C  returns  with  out  confidence  limits  on  function  values 
RETURN 

ASSIGN  414  TO  LAOS 
GO  TO  101 

414  D0  415K=  l.NOB 
TEMP  =  0 

DO  420 1=1, NP 
DO420J=l,NP 
ISUB=K+NOB*(I-l) 

DEBUGl  =  DEL2:(ISUB) 

ISUB  =  K+NOB*(J-l) 

DEBUG2  =  DELZ(ISUB) 

IJ  =  I  +  NP*(J-1) 

DEBUG3  =  D(IJ)/(DIFZ(I)*'raa)*DIFZ(J)*TH(J)) 

420  TEMP  =  TEMP  +  DEBUGl  *  DEBUG2  *  DEBUGS 
TEMP  =  2.0*SQRT(TEMP)*SDEV 
R(K)=F(K)+TEMP 

415  F(K)=F(K)-TEMP 
WRITE(4,  1008) 

IE=0 

DO  425  I=1,NOB,10 
IE=IE-M0 

IF(NOB-IE)  430,435,435 
430  IE=NOB 

435  WRnE(4,  2001)  (R(J),  J  =  I,  E) 

425  WRITE(4,  2006)  (F(J),  J  =  I,  E) 

410  WRITE(4,  1033)  NPROB 
RETURN 

99  WRITE(4, 1034) 

GO  TO  410 

10000FORMAT(38H  NON-LINEAR  ESTIMATION,  PROBLEM  NUMBER  13,// 15, 

1  14H  OBSERVATIONS,  15,  1 IH  PARAMETERS  114,  17H  SCRATCH  REQUIRED) 

1001  FORMAT(/25H  INITIAL  PARAMETER  VALUES  ) 

1002  FORMAT(/54H  PROPORTIONS  USED  IN  CALCULATING  DIFFERENCE 
QUOTENTS  ) 

1003  FORMAT(/25H  INITIAL  SUM  OF  SQUARES  =  E12.4) 

1004  FORMAT(/////45X,  1 3HITERATION  NO.  14) 
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1007  FORMAT(/32H  PARAMETER  VALUES  VIA  REGRESSION  ) 

1008  FORMAT(////54H  APPROXIMATE  CONHDENCE  LIMITS  FOR  EACH  FUNCTION 
VAL 

lUE  ) 

10090FORMAT(/62H  ITERATION  STOPS  -  RELATIVE  CHANGE  IN  EACH  PARAMETER 
LE 

ISSTHAN  E12.4) 

10100FORMAT(/62H  ITERATION  STOPS  -  RELATIVE  CHANGE  IN  SUM  OF  SQUARES 
LE 

ISSTHAN  E12.4) 

101 1  FORMAT(22HO  FUNCTION  VALUES  ) 

1012  FORMAT(////10H  RESIDUALS  ) 

1014  FORMAT(y/24H  VARIANCE  OF  RESIDUALS  =  ,E12.4,1H,I4, 

120H  DEGREES  OF  FREEDOM  ) 

1015  FORMAT(////19H  CORRELATION  MATRIX  ) 

1016  FORMAT(////21H  NORMALIZING  ELEMENTS  ) 

1033  FORMAT(//19H  END  OF  PROBLEM  NO.  13) 

1034  FORMAT(/16H  PARAMETER  ERROR  ) 

10390FORMAT(/71H  INDIVIDUAL  CONFIDENCE  LIMITS  FOR  EACH  PARAMETER  (ON 
LI 

INEAR  HYPOTHESIS)  ) 

10400FORMAT(/9H  LAMBDA  =E10.3,40X,33HSUM  OF  SQUARES  AFTER  REGRESSION 
1E15.7) 

1041  FORMAT(14H  DETERMINANT  =  EI2.4,  6X,  25H  ANGLE  IN  SCALED  COORD.  = 

1  F5.2,  8HDEGREES  ) 

1043  FORMAT(28H  TEST  POINT  SUM  OF  SQUARES  =  El  2.4) 

2001  FORMAT(/10E12.4) 

2006  FORMAT(10E12.4) 

END 

SUBROUTINE  MATIN(A,  NVAR,  B,  NB,  DET) 

DIMENSION  A(NVAR,  1),  B(NVAR,  1) 

PIVOTM  =  A(l,l) 

DET=  1.0 

DO  550 ICOL  =  1,  NVAR 
PIVOT  =  A(ICOL,ICOL) 

PIVOTM  =  AMINl(PIVOT,  PIVOTM) 

DET  =  PIVOT*  DET 

C  DIVIDE  PIVOT  ROW  BY  PIVOT  ELEMENT 
A(ICOL,  ICOL)  =  1.0 
PIVOT=  AMAXl (PIVOT,  l.E-20) 

PIVOT  =  A(ICOL,  ICOL)/PIVOT 
DO  350  L=1,NVAR 
350  A(ICOL,  L)  =  A(ICOL,  L)*PIVOT 
IF(NB  .EQ.  0)  GO  TO  371 
DO  370  L=1,NB 

370  B(ICOL,  L)  =  B(ICOL,  L)*PIVOT 
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C  REDUCE  NON-PIVOT  ROWS 
371  DO550Ll=l,NVAR 

IF(L1  .EQ.  ICOL)  GO  TO  550 
T  =  A(Ll,ICOL) 

A(L1,  ICOL)  =  0. 

DO  450L=1^AR 

450  A(L1,L)  =  A(L1,L)-A(IC0L,L)*T 
IF(NB  .EQ.  0)  GO  TO  550 
DO  500L=1,NB 

500  B(Ll,L)  =  B(Ll,L)-B(ICOL,L)*T 
550  CONTINUE 
RETURN 
END 

SUBROUTINE  GASS60(ITYPE,  NQ,  A,  B,  Q 
DIMENSION  A(NQ),B(NQ),C(NQ,NQ) 

NP  =  NQ 
NR  =  NP/10 
LOW=  1 
LUP=  10 

10  IF( NR)  15,20,30 
15  RETURN 
20  LUP=NP 

IF(LOW  .GT.  LUP)  RETURN 
30  WRITE(4, 500)  (J,J=LOWJ.UP) 

GO  TO  (40,60,80),ITYPE 
40  WRITE(4,  600)  (A(J),J=LOW,LUP) 

GO  TO  100 

60  WRITE(4,  600)  (B(J),J=LOWJLUP) 

GO  TO  40 

80  DO  90  I=LOW,LUP 
90  WRITE(4, 720)  I,(C(J,I),J=LOW,I) 
LOW2=LUP+l 

IF(LOW2  .GT.  NP)  GO  TO  100 
DO  95  I=LOW2,NP 

95  WRITE(4,  720)  I,(C(J,I),J=LOWEUP) 

100  LOW  =  LOW  +  10 
LUP  =  LUP  +  10 
NR  =  NR-  1 
GO  TO  10 

500  FORMAT(/I8,9I12) 

600  FORMAT(10E12.4) 

720  FORMAT(1HO,I3,1X,F7.4,9F12.4) 

1  CONTINUE 
RETURN 
END 

C  THE  MODEL 
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FDSOLVE 

SUBROUTINE  FDSOLVE(TIME,CTIME,ITSTEPS,ANS,XIC, 
.FN,DELTA,CL,STOLD,VS) 

IMPLICIT  NONE 

DOUBLE  PRECISION  ANS(101),XIC(10I), 

.DELTA(10I)W101),K(101),RES(101),A(I01),B(101),TIME,CTIME 
. , RESN, DEEP, RKD, CL, FAST.  ANS2 
INTEGER  ITSTEPS,NFRAC 

DOUBLE  PREQSION  TBCC(30),XQ(30),XCIN(30),FOC,XM,V,PC0 
DOUBLE  PRECISION  X(100),XMR,XKPR,XKR,XFR,DPATH,XSS0,XSR0 
DOUBLE  PRECISION  SNEW(101),SOLD(101),STOLD,ST,VR(30) 
DOUBLE  PRECISION  DELTAT,VS,VA(30),CA(30) 

DIMENSION  IX(100),Y(I00),P(7),DIFF(7),SIGNS(7),SCRAT(1999) 
DIMENSION  IFC(30) 

CHARACTER*8  PNAME(7) 

DIMENSION  VAL(7),IPOINT(7),IVARY(7) 

COMMON /BMOD/FOC,XM,VT>CO,NEFD,TBCC,XQ,XCIN,VR.VA,CA 
COMMON  /BMODA/XMR,XKPR,XKR,XFR,NFRAC,IFLAG 
COMMON  /VALUE/  VAL,IPOINT,IVARY,PNAME,NFIX 
COMMON  /BDAT/  IX,X 
C 

DELTAT=(CTIME-TIME)/ITSTEPS 

C 

C  This  version  of  FDSOLVE  solves  the  compartment 
C  equations  using  backward  Euler  time  differencing 
C  to  avoid  stability  problems.  The  first-order 
C  accuracy  should  be  OK  since  it  matches  the  first-order 
C  spatial  accuracy 
C 

C  Note  a  Gauss-Seidel  solver  is  used  to  solve  the 
C  tiidiagonal  matrix  system  that  results  from  the 
C  Backward  Euler  application 
C 

C  Set  initial  conditions  from  XIC  input 

C  Also  use  this  for  initial  guess  for  Gauss-Seidel  solver 

C 

RKD=VAL(1) 

DEFF=VAL(2) 

DPATH=VAL(3) 

RLAM=VAL(4) 

XSS0=VAL(5) 

FAST=VAL(6) 

DO  10I=1,NFRAC 
SOLD(I)=XIC(I) 

SNEW(I)=XIC(I) 

10  CONTINUE 
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c 

NC0UNT=1 

C 

c 

C  Set  up  constant  used  to  avoid  extra  calculations 
C  see  Dave's  notes 

K(1)=2.0*DEFF*DELTAT/(DELTA(I)**2) 

DO  20  J=2,NFRAC 

K(J)=2.0*DEFF*DELTAT/((DELTA(J-1)+DELTA(J))*DELTA(J)) 
20  CONTINUE 
C 

C  Set  up  matrix  constants 
DO  30J=1,NFRAC-1 
C  Diagonal  term 

A(J)=l.(}t-K(J)+FN(J+l)*K(J+l)/FN(J) 

30  CONTINUE 
C 

A(NFRAC)=  1 .0+K(NFRAC) 

DO  40  J=2,NFRAC 
C  Subdiagonal  term 

B(J)=-K(J)*FN(J)/FN(J-1) 

40  CONTINUE 

C  Note  that  the  super  diagonal  tenn  is  just 
C  equal  to  -k(j),  j=  1  mfrac- 1 
C 


DO  WHILE  (NCOUNT.LE.ITSTEPS) 

C 

C  Gauss-Seidel  iteration 
IGITER=0 
RESN=1.0 

DO  WHILE((RESN.GT.0.0000 1 ). AND.(IGITER.LT.50)) 

C 

SNEW(l)=(K(2)*SNEW(2)+SOLD(l)+K(l)*FN(l) 

*RKD*(1.D0-FAST)*CL)/A(1) 

C 

DO  50  J=2,NFRAC-1 

SNEW(J)=(SOLD(J)+K(J+l)*SNEW(J+l)-B(J)*SNEW(J-l))/A(J) 

C 

50  CONTINUE 

SNEW(NFRAC)=(SOLD(NFRAC)-B(NFRAC)*SNEW(NFRAC-l))/A(NFRAC) 

C 

C  compute  residual 

C 

RESN=0.0 

RES(  1)= A(1  )*SNE  W(  1)-K(2)  *SNE  W(2)- 
SOLD(  1)-K(1)*FN(1)*RKD*CL*(1.D0-FAST) 
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RESN=RESN+(RES(1))**2 

C 

DO  60  J=2,NFRAC-1 
C 

RES(J)=A(J)*SNEW(J)-K(J+1)*SNEW(J+1)+ 

B(J)*SNEW(J-l)-SOLD(J) 

RESN=RESN+(RES(J)**2) 

60  CONTINUE 
C 

RES(NFRAC)=A(NFRAC)*SNEW(NFRAQ+ 
B(NFRAC)*SNEW(NFRAC-l)-SOLD(NFRAC) 
RESN=RESN+(RES(NFRAC)**2) 
RESN=SQRT(RESN) 

IGITER=IGITER+1 

End  Gauss-Seidel  loop 

END  DO 

Add  the  sorbed  concentrations  to  update  cl 

ST=SNEW(1) 

DO  70  J=2NFRAC 
ST=ST+SNEW(J) 

CONTINUE 

Update  the  aqueous  concentration  using  the 
sorbed  concetration  solutions 

CL=CL-i-(XM/(VS+RKD*FAST*XM))*(STOLD-S'i:) 
ANS(NCOUNT)=CL 
ANS2=(ST+FAST*RKD*CL)/CL 
IF  (IFLAG.EQ.0)  ANS(NCOUNT)=ANS2 
C 

DO  80I=1,NFRAC 
SOLD(I)=SNEWa) 

XIC(I)=SNEW(I) 

80  CONTINUE 
C 

STOLD=ST 
NCOUNT=NCOUNT+l 
END  DO 
RETURN 
END 
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SUBROUTINE  INTERP(XT, TIME, CTIME,ANS,ITSTEPS,FR) 
DOUBLE  PRECISION  ANS(101).TIME,CTIME,XT,DELTAT 
DOUBLE  PRECISION  TIMEC 
INTEGER  ITSTEP 
C 

DELTAT=(CnME-TIME)/ITSTEPS 

NCOUNT=0 

TIMEC=TIME 

DO  WHILE  (XT.GT.TIMEC) 

TIMEC=TIMEC+DELTAT 
NCOUNT=NCOUNT+l 
END  DO 

C  FR=(ANS(NCOUNT)+ANS(NCOUNT-1))/2.DO 

FR=(ANS(NCOUNT)-ANS(NCOUNT-l))*(XT-TIMEC+DELTAT) 

./DELTAT+ANS(NCOUNT-l) 

IF  (XT.EQ.TIME)  FR=ANS(1) 

IF  (XT.EQ.CTIME)  FR=ANS(ITSTEPS) 

RETURN 

END 
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INTP 

SUBROUTINE  INTF(RLAM, DELTA, DMAXJ^JFRAC,FN) 

DOUBLE  PRECISION  DELTA(I00),FN(I00),RLAM,DMAX 
INTEGER  NFRAC 
C 

FAC=(RLAM+1)/(DMAX**(RLAM+1)) 

I=I 

X=0.0 

DO  WHILE  a.LT.NFRAC) 

C  Apply  Simpson's  mle 

C  FN(I)=DELTA(I)*FAC/6.0*((DMAX-X)**RLAM+ 

C  .4.0*(DMAX-(X+DELTAa)/2))**RLAM+ 

C  .(DMAX-(X+DELTA(I)))**RLAM) 

FN(I)=(((DMAX-X)**(RLAM+I.DO))-(DMAX-(X+DELTA(I)))** 

.  (RLAM+LD0))/(DMAX**(RLAM+1.D0)) 

X=X+DELTA(I) 

1=1+1 
END  DO 

C  FN(NFRAC)=FAC*DELTA(NFRAO 

FN(NFRAQ=(DELTA(NFRAC)**(RLAM+1.D0))/(DMAX**(RLAM+I.D0)) 

RETURN 

END 
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BATCHI 

PROGRAM  MODEL 

CDLC  SUBROUTINE  MODEL  (NPROB,Pf  ,NOB,NP) 
C  IMPLICIT  NONE 

DOUBLE  PRECISION  XIC(10I),DELTA(101), 

.  FN(101),RKD,CL 
INTEGER  FDOMAIN 


DOUBLE  PRECISION  TBCC(30),XQ(30),XCIN(30),FOC,XM,VPC0 
DOUBLE  PREQSION  X(100),XMR,XKPR.XKR,XFR,DEFF,DPATH,XSSO,XSR0 
DOUBLE  PREaSION  STOLD,VR(30),VS,VA(30),CA(30),Y(100) 

DOUBLE  PRECISION  ANS(101),TIME,CTIME,RLAM,FAST 
DIMENSION  F(100),IFC(30) 

C 

CHARACTER*8  PNAME(7) 

DIMENSION  VAL(7),IPOINT(7),IVARY(7) 

COMMON /BMOD/FOC,XM,V,PCO,NEFD,TBCC,XQ,XCIN,VR.VA,CA 

COMMON  /BMODA/XMR,XKPR,XKR,XFRJSFRAC,IFLAG 

COMMON  A^ALUE/  VAL,IPOINT,IVARY,PNAME,NnX 

COMMON  /BDAT/  IX, X 

CHARACTER*15  INFIL,OUTl,OUT2 

CHARACTER*80  TIT,LrNl  J.IN2,LIN3,LIN4 

NFIX=6 

PNAME(1)='  Kd ' 

PNAME(2)='  Deff ' 

PNAME(3)='  dpath ' 

PNAME(4)='  Lambda ' 

PNAME(5)='  SsO  ’ 

PNAME(6)='  Fast ' 

WRITE(*,*)  'NAME  OF  INPUT  FILE?' 

READ(5,1111)INFIL 
nil  FORMAT(A15) 

WRITE(6,*)NAME  OF  SHORT  OUTPUT  FILE?' 

READ(5,llll)OUTl 

OPEN(UNIT=3,FILE=INHL,STATUS='OLD') 

OPEN(UNIT=10,FILE=OUT1,STATUS='NEW’) 

READ(3,'(A)')  TIT 
WRITE(10,*)TIT 

WRITE(10,*)'  'HME  PREDICTED  CONC  OBSERVED  CONC 
C 

C  Read  parameters  to  be  fitted 
C 

C  RKD  =  partition  coefficient,  KP, 

C  DEFF  =  Coefficient  of  Diffusion,  D_eff  cm''2/sec 
C  DPATH  =  Total  diff.  path  length,  d,  cm 
C  REAM  =  shape  factor  Lambda 
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C  XSSO  =  Initial  dmls  Concentration  in  slow  soil  sites,  Sso,  ug/g 
C  FAST  =  Fraction  of  equilibrium  sorption  sites 
C 

READ(3,'(A)')  LINl 

READ(3,*)  RKD,DEFF,DPATH,RLAM,XSSOf  AST 
C 

C  Read  Constants 
C 

C  FOC  =  Mass  fraction  organic  carbon  on  soil,  foe 
C  XM  =  Mass  of  soil  in  reactor,  g 
C  V  =  Volume  of  solvent  in  reactor  at  time  zero,  mL 
C  PCO  =  initial  concentration  in  liquid,  mg/I 
C  NEFD  =  Number  of  boundary  condition  changes 
C  IFC  =  boundary  condition  number 
C  TBCC  =  time  for  boundary  condition  to  end 
C  XQ  =  Flow  rate,  mLAnin 
C  XCEM  =  Dimensionless  inlet  concentration 
C  XMR  =  Mass  of  reactor,  g 
C  XKPR  =  Linear  partition  coef.  on  reactor,  mL/g 
C  XKR  =  1st  order  mass  transfer  rate  onto  reactor,  min'^-1 
C  XFR  =  Fraction  instaneous  sorbing  reactor  sites 
C  VS(J)  =  volume  of  liquid  in  rx  for  BC  period  j,  ml 
C  VR(J)  =  volume  of  liquid  removed  at  end  of  BC  period!,  ml. 

C  VA(j)  =  volume  of  liquid  added  at  end  of  BC  period  J,  ml. 

C  CA(J)  =  concentration  in  liquid  added  at  end  of  BC  period  Jjmg/l. 
C  XIC  =  sorbed  concentration  in  each  slow  compartment,  ug/g 
C  CL  =  current  liquid  concentration,  mg/1 
C 

READ(3,'(A)')  LIN2 
READ(3,*)  (IVARY(I),I=l,NnX) 

READ(3,'(A)')  LIN3 
READ(3,*)  XM,VJ’C0 
C  READ(3,*)  XMR,XKPR,XKR,XFR 
READ(3,*)  NEFD 
DO  6 1=1, NEFD 

6  READ(3,*)  IFCa),TBCC(I),VR(I),VA(I),CA(I) 

C  Read  numerical  parameter 
C  IFLAG=0  for  sorbed/liquid  concentration 
C  IFLAG=1  for  liquid  concentration 
READ(3,*)NFRAC,IFLAG 
CONTINUE 
C 

C  Read  Data  to  be  fitted 
C 

C  X  =  time,  seconds 
C  Y  =  observation  at  time  X,  see  IFLAG. 
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READ(3,’(A)')  LIN4 

NOB=0 

DO  3  1=1,100 

READ(3,*,END=4)  X(l),Ya) 

3  NOB=NOB+l 

4  CONTINUE 

C 

C  SET  VAL  PARAMETERS  TO  PASS  TO  SUBS 
VAL(1)=RKD 
VAL(3)=DPATH 
VAL(2)=DEFF 
VAL(4)=RLAM 
VAL(5)=XSS0 
VAL(6)=FAST 
C 

CDLC 

C  Begin  new  program  for  batch  solution 
C 
C 
C 

CL=(PCO*V+XM*XSSO*FAST)/(V+RKD*XM*FAST) 

VS=V 

C 

C  Input  the  domain  pardoning.  Delta®  is  the  delta  for 
C  each  compartment. 

C 

DO  101=1,NFRAC 
DELTA(1)=DPATH/NFRAC 
10  CONTINUE 
C 

C  Subroutine  INTF  generates  the  F(i)'s  by  integrating  f(delta) 
C  with  REAM  as  the  shape  parameter 
C 
C 

C  Set  the  initial  condition  within  the  immobile  zone 

C 

C 

CALL  1NTF(RLAM, DELTA, DP  ATH,NFRAC,FN) 

C 

C  STOLD  is  the  sum  of  the  S  for  the  diffusion  compartments 

C 

STOLD=0.0 
DO  20  1=1,NFRAC 
X1C(1)=XSS0*FN(1)*(  1  .DO-FAST) 
STOLD=STOLD+XlC(l) 

20  CONTINUE 
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uuu  uuuuuuu 


Loop  through  each  flow  domain 

FDOMAIN=l 

TIME=0.0 

IF  (IFLAG.EQ.1)  THEN 
WRITE(10,777)TIME,PC0 
ELSE 

WRITE( 1 0,777)TIME,XSS0/PC0 
END  IF 


IT=1 

Note  the  the  final  time  is  X(NOB) 

DO  WHILE  (FDOMAIN.LE.(NEFD+l)) 

ITSTEPS  is  the  number  of  time  steps  per  flow  domain 

ITSTEPS=100 


CnME  is  the  cumulative  solution  time 
CTIME=TBCC(FDOMAIN) 

IF  (FDOMAIN.EQ.(NEFD+l))  CTIME=X(NOB) 

C  Gieck  for  short  flow  domains,  and  adjust  time  steps 
IF  ((CTIME-TIME).LT.10.0)  ITS'EEPS=20 
c 

C  FDSOLVE  computes  the  soln  from  time  to  ctime  for  itsteps  and 
C  passes  the  soln  back  as  array  ANS(ITSTEPS) 
c 

CALL  FDSOLVE(TIME,CnME,ITSTEPS,ANS,XIC, 

.  FN,DELTA,CL,STOLD,VS) 

C 

C  INTERP  interpolates  the  soln  from  FDSOLVE  onto  the  observation 

C  times  and  stores  the  results  in  array  F(nobs) 

C 

C  DO  WHILE  ((X(n:).LE.CTIME).AND.(X(n:).GE.TIME)) 

C  CALL  INTERP(X(IT),TIME,CnME,ANS,ITSTEPS,F(IX)) 

C  IT=IT+1 

C  END  DO 

C  Write  prediction  data  to  output  file 
PnME=TIME 

DELTAT=(CTIME-TIME)/ITSTEPS 
DO  77 11=10, ITSTEPS, 10 
PTIME=PTIME+ 1 0.D0*DELTAT 
WRITE(  10,777)PT1ME,ANS(II) 

777  F0RMAT(1X,2(E16.5)) 

77  CONTINUE 
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non 


c 

TIME=CTIME 

C 

C  Adjust  aqueous  volume  and  new  aqueous  cone. 

C 

C 

c 

CL=(VS*CL+VA(TFDOMAIN)*CA(FDOMAIN)-VR(FDOMAIN)*CL)/ 
.  (VS-VR(FDOMAIN)+VA(FDOMAIN)) 
VS=VS-VR(FDOMAIN)+VA(FDOMAIN) 

FDOMAIN=FDOMAIN+l 


C 

END  DO 


Output  the  results 


1=1 

DO  WHILE  (I.LT.NOB) 
WRnE(10,778)X(I),Ya) 

778  F0RMAT(1X,E16.5,16X,E16.5) 
1=1+1 
END  DO 
STOP 

C  RETURN 
END 
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BATCHMSP 

CDLC  PROGRAM  MODEL 

SUBROUTINE  MODEL  (NPROB,P,F,NOB,NP) 

C  IMPLICIT  NONE 

DOUBLE  PRECISION  XIC(101),DELTA(I01), 

.  FN(101),RKD,CL 
INTEGER  FDOMAIN 
C 

C  Note  for  linking:  I'm  adding  vr(30),va(30),  and  ca(30) 

C  to  common  block  bmod  make  sure  that  this  mode  is  also 
C  made  to  the  caling  program 
C 

DOUBLE  PREQSION  TBCC(30),XQ(30),XCIN(30),FOC,XM,VJPC0 
DOUBLE  PREQSION  X(100),XMR,XKPR,XKR,XFR,DEFF.DPATH,XSS0,XSR0 
DOUBLE  PRECISION  STOLD,VR(30),VS,VA(30),CA(30)TtK2(101) 

DOUBLE  PRECISION  ANS(101),TIME,CTIMETILAM,FAST 
DIMENSION  IX(IOO) 

CDLC 

DIMENSION  F(1)T>(1) 

C 

CHARACTER*8  PNAME(7) 

DIMENSION  VAL(7),IPOINT(7),IVARY(7) 

COMMON /BMOD/FOC,XM,V,PCO,NEFD,TBCC,XQ,XCIN,VR,VA,CA 
COMMON  /BMODAyXMR,XKPR,XKR,XFR,NFRAC,IFLAG 
COMMON  /VALUE/  VAL,IPOINT,IVARY,PNAME JMHX 
COMMON  /BDAT/  IX, X 
C 

C  READ  DATA  FROM  PARAMETERS 
RKD=VAL(1) 

DPATH=VAL(3) 

DEFF=VAL(2) 

RLAM=VAL(4) 

XSS0=VAL(5) 

FAST=VAL(6) 

CDLC  Fake  input  data  to  make  this  a  standalone  program 

CDLC  PC0=0.0 

CDLC  NOB=10 

CDLC  DO100I=lNIOB 

CDLC  X(NOB)=  1.0*1 

CDLC  100  CONTINUE 

CDLC  TBCC(1)=3.0 

CDLC  TBCC(2)=2.0 

CDLC  TBCC(3)=4.0 

CDLC  TBCC(4)=L0 

CDLC  XM=.2 

CDLC  V=.3 
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n  o  o  o  n 


CDLC 

DEFF=1.0 

CDLC 

RKD=1.0 

CDLC 

X(l)=.2 

CDLC 

X(2)=1.2 

CDLC 

X(3)=2.7 

CDLC 

X(4)=3.2 

CDLC 

X(5)=4.8 

CDLC 

X(6)=5.1 

CDLC 

X(7)=8.2 

CDLC 

X(8)=8.9 

CDLC 

X(9)=9.4 

CDLC 

X(10)=10.0 

CDLC 

RLAM=2.0 

CDLC 

C  Begin  new  program  for  batch  solution 
C 

C  NFRAC  is  the  number  of  compartments 
CDLC  NFRAC=100 
C 

C  ED:  RKD  IS  THE  KD  AND  CL  IS  THE  AQUEOUS  CONC 
C  YR  is  the  volume  removed  per  sample,  change 
C  this  later  to  an  array 

C  VS  is  the  current  aqueous  volume,  modified  as 
C  samples  are  taken 

C  VA  is  the  volume  added,  CA  is  the  cone  of  added 
C  liquid. 

C 

C 

CL=(PCO*V+XM*XSSO*FAST)/(V+RKD*XM*FAS'i:) 
VS=V 

Input  the  domain  partioning.  Delta(I)  is  the  delta  for 
each  compartment 

FO  is  the  max  f  value  on  the  f  vs  delta  plot 
F0=(RLAM+  1)/DPATH 
DO  10 1=1, NFRAC 
DELTA(I)=FO/NFRAC 
10  CONTINUE 
C 

C  Subroutine  INTF  generates  the  F(i)'s  by  integrating  delta(l) 

C  with  RLAM  as  the  shape  parameter 

C 

C 

C  Set  the  initial  condition  within  the  immobile  zone 

C 

C 
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CALL  INTF(RLAM,DELTA,DPATH,NFRAC,FN) 

C 

C  STOLD  is  the  sum  of  the  S  for  the  diffusion  compartments 
C  Rk20  is  the  first  order  proportion  constant  for  compartment  j 

STOLD=0.0 
DO  20I=1,NFRAC 
XIC(I)=XSS0*FN(I)*(1.D0-FAST) 
RK2(I)=(DEFF*DELTA(I)**2.D0)/(0.5DO*FN(I)**2.D0) 
STOLD=STOLD+XIC(I) 

20  CONTINUE 
C 

C  Loop  through  each  flow  domain 
C 

FDOMAIN=l 

TIME=0.0 

IT=1 

C  Note  the  the  final  time  is  X(NOB) 

DO  WHILE  (FDOMAIN.LE.(NEFD+l)) 

C 

C  ITSTEPS  is  the  number  of  time  steps  per  flow  domain 
C 

ITS'rEPS=100 

c 

c 

C  CTIME  is  the  cumulative  solution  time 
CTIME=TBCC(FDOMAIN) 

IF(FDOMAIN.EQ.(NEFD+l))  CTIME=X(NOB) 

C  Check  for  short  flow  domains,  and  adjust  time  steps 
IF((CnME-TIME).LT.10.0)  ITSTEPS=20 
c 

C  FDSOLVE  computes  the  soln  fiom  time  to  ctime  for  itsteps  and 

C  passes  the  soln  back  as  array  ANS(ITSTEPS) 
c 

CALL  FDSOLVE(TIME, CTIME, ITSTEPS, ANS,XIC, 

.  FN,DELTA,CL,STOLD,VS,RK2) 

C 

C  ENTERP  interpolates  the  soln  from  FDSOLVE  onto  the  observation 
C  times  and  stores  the  results  in  array  F(nobs) 

C 

DO  WHILE  ((X(ir).LE.CTIME).AND.(X(IT).GE.TlME)) 
CALLINTERP(X(IT),TIME,CTIME,ANS,ITSTEPS,F(IT)) 
IT=IT+1 
END  DO 
TIME=CTIME 
C 

C  Adjust  aqueous  volume  and  new  aqueous  cone. 

C 
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c 

c 

CL=(VS*CL+VA(FDOMAIN)*CA(FDOMAIN)-VR(FDOMAIN)*CL)/ 
.  (VS-VR(FDOMAIN)+VA(FDOMAIN)) 
VS=VS-VR(FDOMAIN)+VA(FDOMAIN) 

FDOMAIN=FDOMAIN+ 1 

C 

END  DO 
C 

C  Output  the  results 
C 

CDLC  1=1 

CDLC  DO  WHILE  (I.LT.NOB) 

CDLC  WRITE(**)LFa) 

CDLC  1=1+1 

CDLC  END  DO 
CDLC  STOP 
RETURN 
END 
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Appendix  D:  MSS  Model  Input  Data 


TEST  FOR  MSSFIT:  Small  Beads  Vial  SBl:  Kp  from  experiment 
RKD  DEEP  DPATH  REAM  XSSO  FAST 
0.425  l.OE-8  0.55  2.0  0.0  0.0 

1  MEANS  FIT  PARAMETER,  0  MEANS  KEEP  IT  FIXED 
010000 

CONSTANTS(Ms,Vli,ai) 

6.272  20.5  0.8910 
9 

1  900.  0.8491  0. 0. 

2  2700.  0.959  0.  0. 

3  5400. 0.878  0.  0. 

4  9000.  0.8602  0.  0. 

5  12600.  0.868  0.  0. 

6  16200.  0.874  0.  0. 

7  21600.  0.84  0.  0. 

8  27000.  0.8591  0.  0. 

9  32400.  0.851  0.  0. 

100  0 

TIME  CONCENTRATION(Cs/Cw) 

900.  0.0343 
2700.  0.0443 
5400.  0.0585 
9000.  0.0789 
12600.  0.0869 
16200.  0.0999 
21600.  0.1265 
27000.  0.1390 
32400.  0.1500 
39600.  0.1644 


86 


TEST  FOR  MSSFIT:  Small  Beads  Vial  SB2:  Kp  from  experiment 
RKD  DEFF  DPATH  REAM  XSSO  FAST 
0.425  l.OE-8  0.55  2.0  0.0  0.0 

1  MEANS  FIT  PARAMETER,  0  MEANS  KEEP  IT  FIXED 
0  1  0000 

CONSTANTS(Ms,Vli,Cli) 

6.272  19.0  0.8910134 
9 

1  900.  0.8166  0.  0. 

2  2700. 0.7362  0.  0. 

3  5400. 0.779  0.  0. 

4  9000.  0.7574  0.  0. 

5  12600.  0.8365  0.  0. 

6  16200.  0.8421  0.  0. 

7  21600.  3.5689  0.  0. 

8  27000. 0.8545  0.  0. 

9  32400. 0.8474  0.  0. 

100  0 

TIME  CONCENTRATION(Cs/Cw) 

900.  0.06808 
2700.  0.0784 
5400.  0.0804 
9000.  0.0962 
12600.  0.1148 
16200.  0.1276 
21600.  0.1505 
27000.  0.2344 
32400.  0.1756 
39600.  0.1946 
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TEST  FOR  MSSFIT;  Small  Beads  Combined,  Kp  from  experiment 
RKD  DEFF  DPATH  REAM  XSSO  FAST 
0.425  l.OE-8  0.55  2.0  0.0  0.0 

1  MEANS  FIT  PARAMETER,  0  MEANS  KEEP  IT  FIXED 
0 10000 

CONSTANTS(Ms,Vli,ai) 

2 

6.272  20.5  0.8910  9.0 
6.272  19.0  0.8910  10.0 

1  900.  0.8491  0.  0. 

2  2700.  0.959  0.  0. 

3  5400.  0.878  0.  0. 

4  9000.  0.8602  0.  0. 

5  12600.  0.868  0.  0. 

6  16200.  0.874  0.  0. 

7  21600.  0.84  0.  0. 

8  27000.  0.8591  0.  0. 

9  32400.  0.851  0.  0. 

10  39600.  0.8358  0.  0. 

1  900.  0.8166  0.  0. 

2  2700.  0.7362  0.  0. 

3  5400.  0.779  0.  0. 

4  9000.  0.7574  0.  0. 

5  12600.  0.8365  0. 0. 

6  16200.  0.8421  0. 0. 

7  21600.  3.5689  0.  0. 

8  27000.  0.8545  0. 0. 

9  32400.  0.8474  0.  0. 

100  0 

TIME  CONCENTRATION(Cs/Cw) 

1  900.  0.0343 
1  2700. 0.0443 
1  5400.  0.0585 
1  9000.  0.0789 
1  12600.  0.0869 
1  16200.  0.0999 
1  21600.  0.1265 
1  27000.  0.1390 
1  32400.  0.1500 
1  39600.  0.1644 

1  90900.  0.1947 

2  900.  0.06808 
2  2700.  0.0784 
2  5400.  0.0804 
2  9000.  0.0962 
2  12600.  0.1148 
2  16200.  0.1276 
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2  21600.  0.1505 
2  27000.  0.2344 
2  32400.  0.1756 
2  39600.  0.1946 


TEST  FOR  MSSFIT:  Large  Beads  Vial  LBl:  Kp  from  experiment 
RKD  DEFF  DPATH  REAM  XSSO  FAST 
0.425  l.OE-8  0.75  2.0  0.0  0.0 

1  MEANS  FIT  PARAMETER,  0  MEANS  KEEP  IT  FIXED 
0  1  0000 

CONSTANTS(Ms,Vli,ai) 

10.603  39.0  0.8027509 
17 

1  660.  0.8644  0.0. 

2  1320.  0.8644  0.  0. 

3  1980.  0.8644  0.  0. 

4  2700.  0.8644  0.  0. 

5  3420.  0.8644  0.  0. 

6  4140.  0.8644  0.  0. 

7  7740.  0.8644  0.  0. 

8  11340.  0.8644  0.  0. 

9  14940.  0.8644  0.  0. 

10  18540.  0.8644  0. 0. 

1122140.  0.86440. 0. 

12  25740.  0.8644  0. 0. 

13  29340.  0.8644  0. 0. 

14  32940.  0.8644  0. 0. 

15  61740.  0.8644  0.  0. 

16  88080.  0.8644  0. 0. 

17  220140.  0.8644  0.  0. 

100  0 

TIME  CONCENTRATION(Cs/Cw) 

660.  0.0596 
1320.  0.0700 
1980.  0.0722 
2700.  0.1147 
3420.  0.1080 
4140.  0.1162 
7740.  0.1523 
11340.  0.1794 
14940.  0.1820 
18540.  0.1873 
22140.  0.2132 
25740.  0.2021 
29340.  0.2305 
32940.  0.2289 
61740.  0.2869 
88080.  0.3241 
220140.  0.4234 
988620.  0.4285 
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TEST  FOR  MSSFIT:  Large  Beads  Vial  LB2:  Kp  from  experiment 
RKD  DEFF  DPATH  REAM  XSSO  FAST 
0.425  lOE-8  0.75  2.0  0.0  0.0 

1  MEANS  FIT  PARAMETER,  0  MEANS  KEEP  IT  FIXED 
010000 

CONSTANTS(Ms,Vli,Cli) 

10.603  39.0  0.8028 
16 

1  660.  0.86440.  0. 

2  1320.  0.8644  0.0. 

3  1980.  0.8644  0. 0. 

4  2640.  0.8644  0. 0. 

5  4140.  0.8644  0. 0. 

6  7740.  0.8644  0. 0. 

711340.  0.8644  0.  0. 

8  14940.  0.8644  0.  0. 

9  18540.  0.8644  0.  0. 

10  22140.0.8644  0.  0. 

11  25740.0.8644  0.0. 

12  29340.  0.8644  0.  0. 

13  32940.  0.8644  0.  0. 

14  61740. 0.8644  0.  0. 

15  88680.0.8644  0.  0. 

16  242880.  0.8644  0.  0. 

100  0 

TIME  CONCENTRATION(Cs/Cw) 

660.  0.0717 
1320.  0.1012 
1980.  0.1269 
2700.  0.1492 
4140.  0.1141 
7740.  0.1540 
11340.  0.1703 
14940.  0.1761 
18540.  0.1762 
22140.  0.1959 
25740.  0.2135 
29340.  0.2161 
32940.  0.2366 
61740.  0.3004 
88680.  0.3247 
242880.  0.4126 
991920.  0.4206 
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TEST  FOR  MSSFIT:  Large  Beads  Combined 

RKD  DEFF  DPATH  REAM  XSSO  FAST 

0.425  l.OE-8  0.75  2.0  0.0  0.0 

1  MEANS  HT  PARAMETER,  0  MEANS  KEEP  IT  FIXED 
010000 

CONSTANTS(Ms,Vli,ai) 

2 

10.603  39.0  0.8028 
17.0 

10.603  39.0  0.8028  16.0 

1  660.  0.8644  0.  0. 

2  1320.  0.8644  0.0. 

3  1980.  0.8644  0.  0. 

4  2700.  0.8644  0.  0. 

5  3420.  0.8644  0. 0. 

6  4140.  0.8644  0.  0. 

7  7740.  0.8644  0.  0. 

8  11340. 0.8644  0.  0. 

9  14940.  0.8644  0.  0. 

10  18540.  0.86440.  0. 

11  22140.  0.86440.0. 

12  25740.  0.8644  0.  0. 

13  29340.  0.86440.  0. 

14  32940.  0.8644  0.0. 

15  61740.  0.8644  0.  0. 

16  88080.  0.8644  0.0. 

17  220140.  0.8644  0. 0. 

1  660.  0.8644  0.  0. 

2  1320.  0.8644  0. 0. 

3  1980.  0.8644  0.  0. 

4  2640.  0.8644  0.  0. 

5  4140.  0.8644  0. 0. 

6  7740.  0.8644  0.  0. 

7  11340.  0.8644  0.  0. 

8  14940.  0.8644  0.  0. 

9  18540.  0.8644  0.  0. 

10  22140.  0.8644  0.  0. 

11  25740.  0.8644  0.  0. 

12  29340.  0.8644  0.  0. 

13  32940.  0.8644  0.  0. 

14  61740.  0.8644  0.  0. 

15  88680.  0.8644  0.  0. 

16  242880.  0.8644  0. 0. 

100  0 

TIME  CONCENTRATION(Cs/Cw) 

1  660.  0.0596 
1  1320.  0.0700 
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1  1980.  0.0722 
1  2700.  0.1147 
1  3420.  0.1080 
1  4140.  0.1162 
1  7740.  0.1523 
1  11340.  0.1794 
1  14940.  0.1820 
1  18540.0.1873 
1  22140.  0.2132 
1  25740.  0.2021 
1  29340.  0.2305 
1  32940.  0.2289 
1  61740.  0.2869 
1  88080. 0.3241 
1  220140.  0.4234 

1  988620.  0.4285 

2  660.  0.0717 
21320.  0.1012 
2  1980.  0.1269 
2  2700.  0.1492 
2  4140.  0.1141 
2  7740.  0.1540 
2  11340.  0.1703 
2  14940.  0.1761 
2  18540.  0.1762 
2  22140.  0.1959 
2  25740.  0.2135 
2  29340.0.2161 
2  32940.  0.2366 
2  61740.0.3004 
2  88680.  0.3247 
2  242880.  0.4126 
2  991920.  0.4206 
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TEST  FOR  MSSFIT;  P  Mixed  Beads  Vial  PMBl:  Kp  from  experiment 
RKD  DEFF  DPATH  REAM  XSSO  FAST 
0.9158  l.OE-8  0.75  2.0  0.0  0.0 

1  MEANS  FIT  PARAMETER,  0  MEANS  KEEP  IT  FIXED 
010000 

CONSTANTS(Ms,Vli,ai) 

13.381  39.0  0.8027509 
18 

1  720.  0.8870  0.  0. 

2  1440.  0.8799  0.  0. 

3  2100.  0.9395  0.  0. 

4  2760.  0.8232  0.  0. 

5  4260.  0.8242  0.  0. 

6  7860. 0.8578  0. 0. 

7  11460.  0.83140.  0. 

8  15060.  0.7947  0.  0. 

9  18660.  0.84290.  0. 

10  22260.  0.8458  0.  0. 

11  27660.  0.8307  0. 0. 

12  33060.  1.0374  0.  0. 

13  72960.  0.8308  0. 0. 

14  111960.  0.8108  0.  0. 

15  167520. 0.9417  0. 0. 

16  530580.  0.8180  0. 0. 

17  1394580. 0.8635  0.  0. 

18  3295500. 0.84600. 0. 

100  0 

TIME  CONCENTRATION(Cs/Cw) 

720.  0.0236 
1440.  0.0082 
2100.  0.0134 
2760.  0.024 
4260.  0.0336 
7860.  0.0257 
11460.  0.0643 
15060.  0.0783 
18660.  0.0619 
22260.  0.0697 
27660.  0.0863 
33060.  0.0952 
72960.  0.166 
111960.  0.2231 
167520.  0.2758 
530580.  0.4789 
1394580.  0.7354 
3295500.  0.9075 
4151460.  0.9158 
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TEST  FOR  MSSFIT:  P  Mixed  Beads  Vial  PMB2:  Kp  from  experiment 
RKD  DEFF  DPATH  REAM  XSSO  FAST 
0.8595  l.E-6  0.75  2.0  0.0  0.0 

1  MEANS  FIT  PARAMETER,  0  MEANS  KEEP  IT  FIXED 
010000 

CONSTANTS(Ms,VU,ai) 

13.381  31.0  0.8019 
18 

1  660.  0.7935  0.  0. 

2  1320.  0.7945  0. 0. 

3  1980.  0.8462  0.  0. 

4  2640.0.9261  0.0. 

5  3900.  0.7889  0. 0. 

6  7500.  0.8163  0.0. 

7  11100.0.8398  0.0. 

8  14700.  0.81210.  0. 

9  18300.  0.8150  0.  0. 

10  21900.  0.8221  0.  0. 

11  27660.  1.1298  0.  0. 

12  33060.  1.0613  0.  0. 

13  69000.0.83800.  0. 

14  108660.  0.8415  0.  0. 

15  151860.  0.8327  0.  0. 

16  514920.  0.8637  0. 0. 

17  1378920. 0.85610.  0. 

18  3279840. 0.9363  0. 0. 

100  0 

TIME  CONCENTRATIONCCs/Cw) 

660.  0.0166 
1320.  0.0184 
1980.  0.0254 
2640.  0.0251 
3900.  0.0333 
7500.  0.0895 
11100.  0.0554 
14700.  0.1077 
18300.  0.0794 
21900.  0.0935 
27660.  0.1419 
33060.  0.1167 
69000.  0.1766 
108660.  0.2594 
151860.  0.3428 
514920.  0.5371 
1378920.  0.7964 
3279840.  0.8671 
4146780.  0.8595 
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TEST  FOR  MSSFIT;  P  Mixed  Beads  Vial  Combined 
RKD  DEFF  DPATH  REAM  XSSO  FAST 
0.88765  l.OE-8  0.75  2.0  0.0  0.0 

1  MEANS  FIT  PARAMETER,  0  MEANS  KEEP  IT  FIXED 
010000 

CONSTANTS(Ms,Vli,ai) 

2 

13.381  39.0  0.8027509  18.0 
13.381  31.0  0.8019  18.0 

1  720.  0.8870  0.  0. 

2  1440.  0.8799  0.  0. 

3  2100.  0.9395  0.  0. 

4  2760.  0.8232  0.  0. 

5  4260.  0.8242  0.  0. 

6  7860.  0.8578  0.  0. 

7  11460.  0.8314  0. 0. 

8  15060.  0.7947  0.  0, 

9  18660.  0.8429  0. 0. 

10  22260.  0.8458  0.  0. 

11  27660.0.8307  0.  0. 

12  33060. 1.03740.  0. 

13  72960.  0.8308  0.  0. 

14  111960.0.8108  0.0. 

15  167520.  0.9417  0.  0. 

16  530580.  0.8180  0.  0. 

17  1394580.  0.8635  0. 0. 

18  3295500. 0.8460  0.  0. 

1  660.  0.7935  0.  0. 

2  1320.  0.7945  0.  0. 

3  1980.  0.8462  0.  0. 

4  2640.  0.9261  0.  0. 

5  3900.  0.7889  0.  0. 

6  7500.  0.8163  0.  0. 

7  11100.  0.8398  0.  0. 

8  14700.  0.8121  0.  0. 

9  18300.  0.8150  0.0. 

10  21900.  0.8221  0.  0. 

11  27660.  1.1298  0.  0. 

12  33060. 1.0613  0.  0. 

13  69000.  0.8380  0.  0. 

14  108660.  0.8415  0.  0. 

15  151860.  0.8327  0.  0. 

16  514920.  0.8637  0.  0. 

17  1378920.  0.8561  0.  0. 

18  3279840.  0.9363  0.  0. 

100  0 

TIME  CONCENTRATION(Cs/Cw) 
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1  720.  0.0236 
1  1440.  0.0082 
1  2100.  0.0134 
1  2760.  0.024 
1  4260.  0.0336 
1  7860.  0.0257 
1  11460.  0.0643 
1  15060.  0.0783 
1  18660.  0.0619 
1  22260.  0.0697 
1  27660.  0.0863 
1  33060.  0.0952 
1  72960.  0.166 
1  111960.0.2231 
1  167520.  0.2758 
1  530580.  0.4789 
1  1394580.  0.7354 
1  3295500.  0.9075 

1  4151460.  0.9158 

2  660.  0.0166 
2  1320.  0.0184 
2  1980.  0.0254 
2  2640.  0.0251 
2  3900.  0.0333 
2  7500.  0.0895 
211100.  0.0554 
2  14700.  0.1077 
2  18300.  0.0794 
2  21900.  0.0935 
2  27660.  0.1419 
2  33060.  0.1167 
2  69000.  0.1766 
2  108660.  0.2594 
2  151860.  0.3428 
2  514920.  0.5371 
2  1378920.  0.7964 
2  3279840.  0.8671 
2  4146780.  0.8595 


TEST  FOR  MSSFIT:  Small  Beads  Vial  SB  1:  Kp  from  experiment 
RKD  DEFF  DPATH  REAM  XSSO  FAST 
0.425  0.12132E-6  0.55  2.0  0.0  0.0 

1  MEANS  FIT  PARAMETER,  0  MEANS  KEEP  IT  FIXED 
001100 

CONSTANTS(Ms,VU,ai) 

6.272  20.5  0.8910 
10 

1  900.  0.8491  0.  0. 

2  2700.  0.959  0.  0. 

3  5400.  0.878  0.  0. 

4  9000.  0.8602  0.  0. 

5  12600.  0.868  0.  0. 

6  16200.  0.874  0.  0. 

7  21600.  0.840.  0. 

8  27000.  0.8591  0.  0. 

9  32400.  0.851  0.  0. 

10  39600. 0.8358  0.  0. 

100  0 

TIME  CONCENTRATION(Cs/Cw) 

900.  0.0343 
2700.  0.0443 
5400.  0.0585 
9000.  0.0789 
12600.  0.0869 
16200.  0.0999 
21600.0.1265 
27000.  0.1390 
32400.0.1500 
39600.  0.1644 
90900.  0.1947 
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TEST  FOR  MSSFIT:  Small  Beads  Vial  SB2:  Kp  from  experiment 
RKD  DEFF  DPATH  REAM  XSSO  FAST 
0.425  0.22469E-06  0.55  2.0  0.0  0.0 

1  MEANS  FIT  PARAMETER,  0  MEANS  KEEP  IT  FIXED 
001100 

CONSTANTS(Ms,Vli,Cli) 

6.272  19.0  0.8910134 
10 

1  900.  0.8166  0.  0. 

2  2700.  0.7362  0.  0. 

3  5400.  0.779  0.  0. 

4  9000.  0.7574  0.  0. 

5  12600.  0.8365  0.  0. 

6  16200.  0.8421  0.  0. 

7  21600.  3.5689  0.  0. 

8  27000.  0.8545  0.  0. 

9  32400.  0.84740.  0. 

10  39600.  0.8471  0.  0. 

100  0 

TIME  CONCENTRATION(Cs/Cw) 

900.  0.06808 
2700.  0.0784 
5400.  0.0804 
9000.  0.0962 
12600.  0.1148 
16200.  0.1276 
21600.  0.1505 
27000.  0.2344 
32400.  0.1756 
39600.  0.1946 
90900.  0.2160 
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TEST  FOR  MSSFIT;  Small  Beads  Combined 

RKD  DEFF  DPATH  RLAM  XSSO  FAST 

0.425  0.16507E-6  0.55  2.0  0.0  0.0 

1  MEANS  FIT  PARAMETER,  0  MEANS  KEEP  IT  FIXED 
001100 

CONSTANTS(Ms,Vli,Cli) 

2 

6.272  20.5  0.8910  10.0 
6.272  19.0  0.8910  10.0 

1  900.  0.8491  0.  0. 

2  2700.  0.959  0.  0. 

3  5400.  0.878  0.  0. 

4  9000.  0.8602  0.  0. 

5  12600.  0.868  0. 0. 

6  16200.  0.874  0.  0. 

7  21600.  0.84  0. 0. 

8  27000.  0.8591  0.  0. 

9  32400.  0.851  0. 0. 

10  39600.  0.8358  0.  0. 

1  900.  0.8166  0.  0. 

2  2700.  0.7362  0.  0. 

3  5400. 0.779  0.  0. 

4  9000.  0.7574  0. 0. 

5  12600. 0.8365  0.  0. 

6  16200. 0.8421  0.  0. 

7  21600.  3.5689  0.  0. 

8  27000.  0.8545  0.  0. 

9  32400.  0.8474  0.  0. 

10  39600.  0.8471  0.  0. 

100  0 

TIME  CONCENTRATION(Cs/Cw) 

1  900.  0.0343 
1  2700.  0.0443 
1  5400.  0.0585 
1  9000.  0.0789 
1  12600.  0.0869 
1  16200.  0.0999 
1  21600.  0.1265 
1  27000.  0.1390 
1  32400.  0.1500 
1  39600.  0.1644 

1  90900.  0.1947 

2  900.  0.06808 
2  2700.  0.0784 
2  5400.  0.0804 
2  9000.  0.0962 
2  12600.0.1148 
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2  16200.  0.1276 
2  21600.  0.1505 
2  27000.  0.2344 
2  32400.  0.1756 
2  39600.  0.1946 
2  90900.  0.2160 


TEST  FOR  MSSFIT:  Large  Beads  Vial  LB  1 :  Kp  from  experiment 
RKD  DEFF  DPATH  RLAM  XSSO  FAST 
0.425  0.89949E-6  0.75  2.0  0.0  0.0 

1  MEANS  RT  PARAMETER,  0  MEANS  KEEP  IT  FIXED 
001100 

CONSTANTS(Ms,Vli,Cli) 

10.603  39.0  0.8027509 
17 

1  660.  0.8644  0.  0. 

2  1320. 0.8644  0.  0. 

3  1980.  0.8644  0.  0. 

4  2700.  0.8644  0.  0. 

5  3420. 0.8644  0.  0. 

6  4140.  0.8644  0.  0. 

7  7740.  0.8644  0.  0. 

8  11340.  0.86440.  0. 

9  14940.  0.8644  0.  0. 

10  18540.  0.8644  0.  0. 

11  22140.  0.8644  0.  0. 

12  25740.  0.8644  0. 0. 

13  29340.  0.8644  0.  0. 

14  32940.  0.8644  0.  0. 

15  61740.  0.8644  0. 0. 

16  88080.  0.8644  0.  0. 

17  220140.  0.8644  0.  0. 

100  0 

TIME  CONCENTRATION(Cs/Cw) 

660.  0.0596 
1320.  0.0700 
1980.  0.0722 
2700.  0.1147 
3420.  0.1080 
4140.  0.1162 
7740.  0.1523 
11340.  0.1794 
14940.  0.1820 
18540.  0.1873 
22140.  0.2132 
25740.  0.2021 
29340.  0.2305 
32940.  0.2289 
61740.  0.2869 
88080.  0.3241 
220140.  0.4234 
988620.  0.4285 
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TEST  FOR  MSSFIT:  Large  Beads  Vial  LB2:  Kp  from  experiment 
RKD  DEFF  DPATH  REAM  XSSO  FAST 
0.425  0.92897E-6  0.75  2.0  0.0  0.0 

1  MEANS  HT  PARAMETER,  0  MEANS  KEEP  IT  FIXED 
001100 

CONSTANTS(Ms,Vli,ai) 

10.603  39.0  0.8028 
16 

1  660.  0.8644  0.  0. 

2  1320.  0.8644  0.  0. 

3  1980.  0.8644  0.  0. 

4  2640.  0.8644  0.  0. 

5  4140.  0.8644  0.  0. 

6  7740.  0.8644  0.  0. 

7  11340.  0.8644  0. 0. 

8  14940.  0.8644  0.0. 

9  18540.  0.8644  0.0. 

10  22140.  0.8644  0. 0. 

1125740.0.8644  0.0. 

12  29340.  0.8644  0. 0. 

13  32940.  0.8644  0. 0. 

14  61740.  0.8644  0. 0. 

15  88680.  0.8644  0.  0. 

16  242880.  0.8644  0.0. 

100  0 

TIME  CONCENTRATION(Cs/Cw) 

660.  0.0717 
1320.  0.1012 
1980.  0.1269 
2700.  0.1492 
4140.  0.1141 
7740.  0.1540 
11340.  0.1703 
14940.  0.1761 
18540.  0.1762 
22140.  0.1959 
25740.  0.2135 
29340.  0.2161 
32940.  0.2366 
61740.  0.3004 
88680.  0.3247 
242880.  0.4126 
991920.  0.4206 
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TEST  FOR  MSSFIT:  Large  Beads  Combined 

RKD  DEFF  DPATH  REAM  XSSO  FAST 

0.425  0.91935E-6  0.75  2.0  0.0  0.0 

1  MEANS  FIT  PARAMETER,  0  MEANS  KEEP  IT  FIXED 
001100 

CONSTANTS(Ms.VH,CK) 

2  10.603  39.0  0.8028 

17.0  10.603  39.0  0.8028  16.0 

1  660.  0.8644  0.  0. 

2  1320.  0.8644  0. 0. 

3  1980.  0.8644  0.  0. 

4  2700.  0.8644  0. 0. 

5  3420.  0.8644  0.  0. 

6  4140.  0.8644  0. 0. 

7  7740.  0.8644  0. 0. 

8  11340.0.8644  0.  0. 

9  14940.0.8644  0.  0. 

10  18540.  0.8644  0.  0. 

11  22140.0.8644  0.  0. 

12  25740.  0.8644  0.  0. 

13  29340.  0.8644  0.  0. 

14  32940.0.8644  0.  0. 

15  61740.  0.8644  0.  0. 

16  88080.  0.86440.  0. 

17  220140.  0.8644  0.  0. 

1  660.  0.8644  0. 0. 

2  1320.  0.8644  0.0. 

3  1980.  0.8644  0. 0. 

4  2640.  0.8644  0.  0. 

5  4140.0.8644  0. 0. 

6  7740.  0.8644  0.  0. 

7  11340.0.8644  0.  0. 

8  14940.0.8644  0.  0. 

9  18540.0.8644  0.  0. 

10  22140.0.8644  0.  0. 

11  25740.  0.8644  0.  0. 

12  29340.0.8644  0.  0. 

13  32940.  0.8644  0.  0. 

14  61740.0.8644  0.  0. 

15  88680.  0.86440.  0. 

16  242880.  0.8644  0.  0. 

100  0 

TIME  CONCENTRATION(Cs/Cw) 

1  660.  0.0596 
1  1320.  0.0700 
1  1980.  0.0722 
1  2700.  0.1147 
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1  3420.  0.1080 
14140.  0.1162 
1  7740.  0.1523 
1  11340.  0.1794 
1  14940.  0.1820 
1  18540.  0.1873 
1  22140.  0.2132 
1  25740.  0.2021 
1  29340.  0.2305 
1  32940.  0.2289 
1  61740.  0.2869 
1  88080.  0.3241 
1  220140.  0.4234 

1  988620.  0.4285 

2  660.  0.0717 
2  1320.  0.1012 
2  1980.  0.1269 
2  2700.  0.1492 
2  4140.  0.1141 
2  7740.  0.1540 
2  11340.  0.1703 
2  14940.  0.1761 
2  18540.  0.1762 
2  22140.  0.1959 
2  25740.  0.2135 
2  29340.  0.2161 
2  32940.  0.2366 
2  61740.  0.3004 
2  88680.  0.3247 
2  242880.  0.4126 
2  991920.  0.4206 


TEST  FOR  MSSFIT:  P  Mixed  Beads  Vial  PMB 1;  Kp  from  experiment 
RKD  DEFF  DPATH  REAM  XSSO  FAST 
0.888  1.16E-7  5.5  3.5  0.0  0.0 

1  MEANS  FIT  PARAMETER,  0  MEANS  KEEP  IT  FIXED 
00  1100 

CONSTANTS(Ms,Vli,ai) 

13.381  39.0  0.8027509 
18 

1  720.  0.8870  0.  0. 

2  1440.  0.8799  0.  0. 

3  2100.  0.9395  0.  0. 

4  2760.  0.8232  0.  0. 

5  4260.  0.82420.  0. 

6  7860.  0.8578  0.  0. 

7  11460.  0.83140. 0. 

8  15060.  0.7947  0. 0. 

9  18660.  0.8429  0.  0. 

10  22260. 0.8458  0.  0. 

1127660.  0.8307  0.  0. 

12  33060. 1.0374  0.  0. 

13  72960.  0.8308  0.  0. 

14111960.  0.8108  0.  0. 

15  167520.  0.9417  0.  0. 

16  530580. 0.8180  0.  0. 

17  1394580.  0.8635  0. 0. 

18  3295500.  0.84600.  0. 

100  0 

TIME  CONCENTRATION(Cs/Cw) 

720.  0.0236 
1440.  0.0082 
2100.  0.0134 
2760.  0.024 
4260.  0.0336 
7860.  0.0257 
11460.  0.0643 
15060.  0.0783 
18660.  0.0619 
22260.  0.0697 
27660.  0.0863 
33060.  0.0952 
72960.  0.166 
111960.  0.2231 
167520.  0.2758 
530580.  0.4789 
1394580.  0.7354 
3295500.  0.9075 
4151460.  0.9158 
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TEST  FOR  MSSFIT:  P  Mixed  Beads  Vial  PMB2:  Kp  from  experiment 
RKD  DEFF  DPATH  RLAM  XSSO  FAST 
0.888  1.16E-7  5.0  3.5  0.0  0.0 

1  MEANS  HT  PARAMETER,  0  MEANS  KEEP  IT  FIXED 
001  100 

CONSTANTS(Ms,VH,Cli) 

13.381  31.0  0.8019 
18 

1  660.  0.7935  0.  0. 

2  1320.  0.7945  0.  0. 

3  1980.  0.8462  0.  0. 

4  2640.  0.9261  0.  0. 

5  3900.  0.7889  0.  0. 

6  7500.  0.8163  0.  0. 

7  11100. 0.8398  0. 0. 

8  14700.0.8121  0.0. 

9  18300.  0.8150  0. 0. 

10  21900.  0.8221  0.  0. 

11  27660. 1.1298  0.  0. 

12  33060. 1.0613  0.  0. 

13  69000. 0.8380  0.  0. 

14  108660.  0.8415  0.  0. 

15  151860.  0.83270.0. 

16  514920. 0.8637  0.0. 

17  1378920.  0.8561  0. 0. 

18  3279840.  0.9363  0. 0. 

100  0 

TIME  CONCENTRATION(Cs/Cw) 

660.  0.0166 
1320.  0.0184 
1980.  0.0254 
2640.  0.0251 
3900.  0.0333 
7500.  0.0895 
11100.  0.0554 
14700.  0.1077 
18300.  0.0794 
21900.  0.0935 
27660.  0.1419 
33060.  0.1167 
69000.  0.1766 
108660.  0.2594 
151860.  0.3428 
514920.  0.5371 
1378920.  0.7964 
3279840.  0.8671 
4146780.  0.8595 
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TEST  FOR  MSSFIT:  P  Mixed  Beads  Vial  Combined 
RKD  DEFF  DPATH  REAM  XSSO  FAST 
0.888  1.16E-7  5.0  3.5  0.0  0.0 

1  MEANS  HT  PARAMETER,  0  MEANS  KEEP  IT  FIXED 
001100 

CONSTANTS(Ms,V]i,Cli) 

2 

13.381  39.0  0.8027509  18.0 
13.381  31.0  0.8019  18.0 

1  720.  0.8870  0.  0. 

2  1440.  0.8799  0. 0. 

3  2100.  0.9395  0. 0. 

4  2760.  0.8232  0. 0. 

5  4260.  0.8242  0. 0. 

6  7860. 0.8578  0. 0. 

7  11460.0.83140.  0. 

8  15060. 0.7947  0.  0. 

9  18660. 0.8429  0.  0. 

10  22260.  0.8458  0.  0. 

11  27660.  0.8307  0.  0. 

12  33060. 1.0374  0.  0. 

13  72960.  0.8308  0.  0. 

14  111960.0.8108  0.0. 

15  167520.  0.9417  0. 0. 

16  530580.  0.8180  0. 0. 

17  1394580. 0.8635  0.  0. 

18  3295500. 0.8460  0.  0. 

1  660.  0.7935  0. 0. 

2  1320.  0.7945  0. 0. 

3  1980.  0.8462  0. 0. 

4  2640.  0.9261  0. 0. 

5  3900.  0.7889  0. 0. 

6  7500.  0.8163  0.0. 

7  11100.0.8398  0.  0. 

8  14700.0.8121  0.  0. 

9  18300.0.8150  0.  0. 

10  21900.  0.8221  0.  0. 

11  27660.  1.1298  0.  0. 

12  33060.  1.0613  0.  0. 

13  69000.  0.8380  0.  0. 

14  108660.  0.8415  0. 0. 

15  151860.  0.8327  0. 0. 

16  514920.  0.8637  0.0. 

17  1378920. 0.8561  0.  0. 

18  3279840.  0.9363  0.  0. 

100  0 

TIME  CONCENTRATION(Cs/Cw) 
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1  720.  0.0236 
1  1440.  0.0082 
1  2100.  0.0134 
1  2760.  0.024 
1  4260.  0.0336 
1  7860.  0.0257 
1  11460.  0.0643 
1  15060.  0.0783 
1  18660.  0.0619 
1  22260.  0.0697 
1  27660.  0.0863 
1  33060.  0.0952 
1  72960.  0.166 
1  111960.  0.2231 
1  167520.  0.2758 
1  530580.  0.4789 
1  1394580.  0.7354 
1  3295500.  0.9075 

1  4151460.  0.9158 

2  660.  0.0166 
21320.  0.0184 
2  1980.  0.0254 
2  2640.  0.0251 
2  3900. 0.0333 
2  7500.  0.0895 
211100.  0.0554 
2  14700.  0.1077 
2  18300.  0.0794 
2  21900.  0.0935 
2  27660.  0.1419 
2  33060.  0.1167 
2  69000.  0.1766 
2  108660.  0.2594 
2  151860.  0.3428 
2  514920.  0.5371 
2  1378920.  0.7964 
2  3279840.  0.8671 
2  4146780.  0.8595 
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